Analysis: Microhabitats shape diversity-productivity relationships in freshwater bacterial communities

Load Libraries

library(ggplot2)  
library(devtools)
library(phyloseq)
library(picante)
library(tidyr)
library(forcats) 
library(dplyr)
library(tibble)
library(cowplot)
library(picante)    # Will also include ape and vegan 
library(car)        # For residual analysis
library(sandwich)   # for vcovHC function in post-hoc test
library(MASS)       # For studres in plot_residuals function
library(caret)      # For cross validation
library(pander)     # For pretty tables 
library(glmnet)     # For lasso and ridge regressions
library(broom)      # pvalues out of multiple regression
library(DT)         # fancy HTML table output
source("code/Muskegon_functions.R")
source("code/set_colors.R")  

Load data

# Loads a phyloseq object named otu_merged_musk_pruned)
load("data/surface_PAFL_otu_pruned_raw.RData")
# The name of the phyloseq object is: 
surface_PAFL_otu_pruned_raw 

## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 7806 taxa and 24 samples ]
## sample_data() Sample Data:       [ 24 samples by 44 sample variables ]
## tax_table()   Taxonomy Table:    [ 7806 taxa by 8 taxonomic ranks ]
## phy_tree()    Phylogenetic Tree: [ 7806 tips and 7804 internal nodes ]

# Remove doubletons!
surface_PAFL_otu_pruned_rm2 <- prune_taxa(taxa_sums(surface_PAFL_otu_pruned_raw) > 2, surface_PAFL_otu_pruned_raw) 
surface_PAFL_otu_pruned_rm2

## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 2979 taxa and 24 samples ]
## sample_data() Sample Data:       [ 24 samples by 44 sample variables ]
## tax_table()   Taxonomy Table:    [ 2979 taxa by 8 taxonomic ranks ]
## phy_tree()    Phylogenetic Tree: [ 2979 tips and 2977 internal nodes ]

# Remove tree for computational efficiency 
surface_PAFL_otu_pruned_notree_rm2 <- phyloseq(tax_table(surface_PAFL_otu_pruned_rm2), otu_table(surface_PAFL_otu_pruned_rm2), sample_data(surface_PAFL_otu_pruned_rm2)) 

# Gather the metadata in a dataframe to play with 
metadata <- data.frame(sample_data(surface_PAFL_otu_pruned_notree_rm2)) %>%
      mutate(fraction = factor(fraction, levels = c("WholePart","WholeFree")),
         lakesite = factor(lakesite,  levels = c("MOT", "MDP", "MBR", "MIN")),
         fraction = fct_recode(fraction, Particle = "WholePart", Free = "WholeFree"),
         lakesite = fct_recode(lakesite, Outlet = "MOT", Deep = "MDP", Bear = "MBR", River = "MIN"))
row.names(metadata) <- metadata$norep_filter_name

# Replace the sample data 
sample_data(surface_PAFL_otu_pruned_notree_rm2) <- metadata

# Create a dataframe for environmental data only
environmental_data <- metadata %>%
  dplyr::select(Temp_C:DO_percent, -BGA_cellspermL, -SRP_ugperL, fraction, norep_filter_name, -DO_percent) %>%
  dplyr::filter(fraction == "Free") %>%
  dplyr::select(-fraction) %>%
  tibble::column_to_rownames(var = "norep_filter_name") %>%
  dplyr::rename(Temp = Temp_C, SPC = SpCond_uSpercm,
         TDS = TDS_mgperL, ORP = ORP_mV,
         Chla = Chl_Lab_ugperL, Cl = Cl_mgperL,
         SO4 = SO4_mgperL, NO3 = NO3_mgperL,
         NH3 = NH3_mgperL, TKN = TKN_mgperL,
         TP = TP_ugperL, Alk = Alk_mgperL,
         DO = DO_mgperL) %>%
  as.matrix()

#### PCA DATA 
# Scale the data so their variances are the same!
scaled_enviro <- scale(environmental_data)

# Sanity Checks: check that we get mean of 0 and sd of 1
apply(scaled_enviro, 2, mean)

##          Temp           SPC           TDS            pH           ORP          Chla            Cl           SO4           NO3           NH3           TKN            TP           Alk            DO 
## -3.151262e-16 -1.170392e-15 -8.881784e-16 -2.063155e-15  4.614635e-18  2.081668e-17 -2.937375e-16  3.932311e-16  1.237526e-16 -1.155579e-17  7.400583e-17 -5.090894e-17  1.454855e-15  7.401939e-17

apply(scaled_enviro, 2, sd)

## Temp  SPC  TDS   pH  ORP Chla   Cl  SO4  NO3  NH3  TKN   TP  Alk   DO 
##    1    1    1    1    1    1    1    1    1    1    1    1    1    1

# Run a principal component analysis via the vegan package
pca_environ <- rda(scaled_enviro) 

# Pull out all of the PCA scores into a dataframe
pca_scores_df <- summary(pca_environ)$sites %>%
  as.data.frame() %>%
  tibble::rownames_to_column(var = "norep_filter_name") %>%
  mutate(norep_water_name = paste(substr(norep_filter_name, 1, 4), substr(norep_filter_name, 6, 8), sep = "")) %>%
  dplyr::select(-norep_filter_name)

# Combine the above dataframe with the rest of the metadata
metadata_pca <- metadata %>%
  mutate(norep_water_name = paste(substr(norep_filter_name, 1, 4), substr(norep_filter_name, 6, 8), sep = "")) %>%
  left_join(pca_scores_df, by = "norep_water_name")

Removing the OTUs that are doubletons across the entire dataset leaves a total of 2979 OTUs, a decrease from 7806. This will help decrease the computational load for phylogenetic analysis and also further remove potential sequencing artifacts.

Calculate Community Diversity

# Set the seed for reproducibility
set.seed(777)

# Calculate the alpha diversity with 100 repsampling events
alpha_calc <- calc_alpha_diversity(physeq = surface_PAFL_otu_pruned_notree_rm2)

# What was the minimum sample size? 
min(sample_sums(surface_PAFL_otu_pruned_notree_rm2)) - 1

## [1] 6489

# Put it altogether in a dataframe 
otu_alphadiv <- calc_mean_alphadiv(physeq = surface_PAFL_otu_pruned_notree_rm2,
                           richness_df = alpha_calc$Richness, 
                           evenness_df = alpha_calc$Inverse_Simpson, 
                           shannon_df = alpha_calc$Shannon) %>%
  mutate(fraction = factor(fraction, levels = c("Particle","Free")),
         lakesite = factor(lakesite,  levels = c("Outlet", "Deep", "Bear", "River")),
         measure = factor(measure, levels = c("Richness",  "Shannon_Entropy", "Inverse_Simpson", "Simpsons_Evenness")),
         norep_water_name = paste(substr(norep_filter_name, 1, 4), substr(norep_filter_name, 6, 8), sep = ""))

otu_alphadiv <- otu_alphadiv %>%
  left_join(pca_scores_df, by = "norep_water_name")

# Make subsetted dataframes for each diversity metric
richness <- filter(otu_alphadiv, measure == "Richness")
shannon <- filter(otu_alphadiv, measure == "Shannon_Entropy")
invsimps <- filter(otu_alphadiv, measure == "Inverse_Simpson")
simpseven <- filter(otu_alphadiv, measure == "Simpsons_Evenness")

Calculate the Phylogenetic Diversity

# Read in the tree
RAREFIED_rm2_fasttree <- read.tree(file = "data/PhyloTree/newick_tree_rm2_rmN.tre")
  
# Load in data that has doubletons removed 
load("data/PhyloTree/surface_PAFL_otu_pruned_RAREFIED_rm2.RData")
surface_PAFL_otu_pruned_RAREFIED_rm2

## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 1891 taxa and 24 samples ]
## sample_data() Sample Data:       [ 24 samples by 70 sample variables ]
## tax_table()   Taxonomy Table:    [ 1891 taxa by 8 taxonomic ranks ]
## phy_tree()    Phylogenetic Tree: [ 1891 tips and 1889 internal nodes ]

# Create the OTU table for picante 
surface_PAFL_RAREFIED_rm2_otu <- matrix(otu_table(surface_PAFL_otu_pruned_RAREFIED_rm2), nrow = nrow(otu_table(surface_PAFL_otu_pruned_RAREFIED_rm2)))
rownames(surface_PAFL_RAREFIED_rm2_otu) <- sample_names(surface_PAFL_otu_pruned_RAREFIED_rm2)
colnames(surface_PAFL_RAREFIED_rm2_otu) <- taxa_names(surface_PAFL_otu_pruned_RAREFIED_rm2)
    
  
## Calculate input for SES_MPD  
# Convert the abundance data to standardized abundanced vegan function `decostand' , NOTE: method = "total"
otu_decostand_total <- decostand(surface_PAFL_RAREFIED_rm2_otu, method = "total")
# check total abundance in each sample
apply(otu_decostand_total, 1, sum)

## MBREJ515 MBREJ715 MBREJ915 MBREK515 MBREK715 MBREK915 MDPEJ515 MDPEJ715 MDPEJ915 MDPEK515 MDPEK715 MDPEK915 MINEJ515 MINEJ715 MINEJ915 MINEK515 MINEK715 MINEK915 MOTEJ515 MOTEJ715 MOTEJ915 MOTEK515 MOTEK715 MOTEK915 
##        1        1        1        1        1        1        1        1        1        1        1        1        1        1        1        1        1        1        1        1        1        1        1        1

# check for mismatches/missing species between community data and phylo tree
RAREFIED_rm2_matches <- match.phylo.comm(RAREFIED_rm2_fasttree, otu_decostand_total)
# the resulting object is a list with $phy and $comm elements.  replace our
# original data with the sorted/matched data
phy_RAREFIED_rm2 <- RAREFIED_rm2_matches$phy
comm_RAREFIED_rm2 <- RAREFIED_rm2_matches$comm

# Calculate the phylogenetic distances
phy_dist_RAREFIED_rm2 <- cophenetic(phy_RAREFIED_rm2)

  
## Calculate SES_MPD
###################################### INDEPENDENT SWAP ############################################
# calculate standardized effect size mean pairwise distance (ses.mpd)
unweighted_sesMPD_indepswap_RAREFIED_rm2 <- ses.mpd(comm_RAREFIED_rm2, phy_dist_RAREFIED_rm2, null.model = "independentswap", 
                                     abundance.weighted = FALSE, runs = 999)

WEIGHTED_sesMPD_indepswap_RAREFIED_rm2 <- ses.mpd(comm_RAREFIED_rm2, phy_dist_RAREFIED_rm2, null.model = "independentswap", 
                                     abundance.weighted = TRUE, runs = 999)


# Gather div info
rich_df <- filter(otu_alphadiv, measure == "Richness") %>%
  dplyr::select(norep_filter_name, mean, sd, frac_bacprod, SD_frac_bacprod, fracprod_per_cell_noinf)

invsimps_df <- filter(otu_alphadiv, measure == "Inverse_Simpson") %>%
  dplyr::select(norep_filter_name, mean, sd, frac_bacprod, SD_frac_bacprod, fracprod_per_cell_noinf)


# Prepare to be merged with each other 
unweighted_df <- unweighted_sesMPD_indepswap_RAREFIED_rm2 %>%
  rownames_to_column("names") %>%
  mutate(phylo_measure = "Unweighted_SESMPD") %>%
  make_metadata_norep() %>%
  mutate(fraction = fct_recode(fraction, Particle = "WholePart", Free = "WholeFree")) %>%
  rename(norep_filter_name = names) %>%
  left_join(rich_df, by = "norep_filter_name") %>%
  mutate(lakesite = factor(lakesite, levels = c("MOT", "MDP","MBR", "MIN")),
         lakesite = fct_recode(lakesite, Outlet = "MOT", Deep = "MDP", Bear = "MBR", River = "MIN"))
  
weighted_df <- WEIGHTED_sesMPD_indepswap_RAREFIED_rm2 %>%
  rownames_to_column("names") %>%
  mutate(phylo_measure = "Weighted_SESMPD") %>%
  make_metadata_norep() %>%
  mutate(fraction = fct_recode(fraction, Particle = "WholePart", Free = "WholeFree")) %>%
  rename(norep_filter_name = names) %>%
  left_join(invsimps_df, by = "norep_filter_name") %>%
  mutate(lakesite = factor(lakesite, levels = c("MOT", "MDP","MBR", "MIN")),
         lakesite = fct_recode(lakesite, Outlet = "MOT", Deep = "MDP", Bear = "MBR", River = "MIN"))

Main Figures

Have a legend for all plots

legend_plot <- ggplot(metadata, aes(y = frac_bacprod, x = fraction, fill = fraction, shape = season)) + 
  geom_jitter(size = 3, width = 0.2) + 
  scale_fill_manual(values = fraction_colors, guide = guide_legend(override.aes = list(shape = 22, size = 4))) +
  scale_shape_manual(values = season_shapes) +
  theme(legend.position = "bottom", axis.title.x = element_blank(),
        legend.title = element_blank())

season_legend <- get_legend(legend_plot)

Figure 1

######################################################### Fraction ABUNDANCe 
frac_abund_wilcox <- wilcox.test(log10(as.numeric(fraction_bac_abund)) ~ fraction, data = metadata)
frac_abund_wilcox

## 
##  Wilcoxon rank sum test
## 
## data:  log10(as.numeric(fraction_bac_abund)) by fraction
## W = 0, p-value = 1.479e-06
## alternative hypothesis: true location shift is not equal to 0

metadata %>%
  group_by(fraction) %>%
  summarize(mean(as.numeric(fraction_bac_abund)))

## # A tibble: 2 x 2
##   fraction `mean(as.numeric(fraction_bac_abund))`
##     <fctr>                                  <dbl>
## 1 Particle                               41168.88
## 2     Free                              734522.25

# Make a data frame to draw significance line in boxplot (visually calculated)
dat1 <- data.frame(a = c(1.1,1.1,1.9,1.9), b = c(6.45,6.5,6.5,6.45)) # WholePart vs WholeFree

poster_a <- 
  ggplot(filter(metadata, norep_filter_name != "MOTEJ515"), 
       aes(y = log10(as.numeric(fraction_bac_abund)), x = fraction)) +
  geom_jitter(size = 3, aes(fill = fraction, shape = season), width = 0.2) + 
  geom_boxplot(alpha = 0.5, outlier.shape = NA, aes( fill = fraction)) +
  ylab("Log10(Bacterial Counts) \n (cells/mL)") +
  scale_fill_manual(values = fraction_colors) +
  scale_shape_manual(values = season_shapes) +
  ##### Particle vs free cell abundances 
  geom_path(data = dat1, aes(x = a, y = b), linetype = 1, color = "gray40") +
  annotate("text", x=1.5, y=6.55, size = 8, color = "gray40", label= paste("***")) +
  annotate("text", x=1.5, y=6.4, size = 3.5, color = "gray40",
           label= paste("p =", round(frac_abund_wilcox$p.value, digits = 6))) +
  theme(legend.position = "bottom",
        legend.title = element_blank(),
        axis.title.x = element_blank())



######################################################### TOTAL PRODUCTION 
totprod_wilcox <- wilcox.test(frac_bacprod ~ fraction, data = metadata)
totprod_wilcox

## 
##  Wilcoxon rank sum test
## 
## data:  frac_bacprod by fraction
## W = 33, p-value = 0.02418
## alternative hypothesis: true location shift is not equal to 0

metadata %>%
  group_by(fraction) %>%
  summarize(mean(frac_bacprod))

## # A tibble: 2 x 2
##   fraction `mean(frac_bacprod)`
##     <fctr>                <dbl>
## 1 Particle             9.958333
## 2     Free            24.058333

# Make a data frame to draw significance line in boxplot (visually calculated)
dat2 <- data.frame(a = c(1.1,1.1,1.9,1.9), b = c(71,72,72,71)) # WholePart vs WholeFree


poster_b <- ggplot(metadata, aes(y = frac_bacprod, x = fraction)) + 
  geom_jitter(size = 3, aes(fill = fraction, shape = season), width = 0.2) + 
  geom_boxplot(alpha = 0.5, outlier.shape = NA, aes( fill = fraction)) +
  scale_fill_manual(values = fraction_colors, guide = FALSE) +
  scale_shape_manual(values = season_shapes) +
  scale_y_continuous(limits = c(0, 78), expand = c(0,0), breaks = seq(0, 75, by = 15)) +
  ylab("Community Production \n (μgC/L/day)") +
  ##### Particle vs free bulk production  
  geom_path(data = dat2, aes(x = a, y = b), linetype = 1, color = "gray40") +
  annotate("text", x=1.5, y=73, size = 8, color = "gray40", label= paste("*")) +
  annotate("text", x=1.5, y=69, size = 3.5, color = "gray40",
           label= paste("p =", round(totprod_wilcox$p.value, digits = 3))) +
  theme(legend.position = "none",
        axis.title.x = element_blank())


######################################################### TOTAL PRODUCTION 
percellprod_wilcox <- wilcox.test(log10(fracprod_per_cell) ~ fraction, data = filter(metadata, norep_filter_name != "MOTEJ515"))
percellprod_wilcox

## 
##  Wilcoxon rank sum test
## 
## data:  log10(fracprod_per_cell) by fraction
## W = 125, p-value = 6.656e-05
## alternative hypothesis: true location shift is not equal to 0

filter(metadata, norep_filter_name != "MOTEJ515") %>%
  group_by(fraction) %>%
  summarize(mean(fracprod_per_cell))

## # A tibble: 2 x 2
##   fraction `mean(fracprod_per_cell)`
##     <fctr>                     <dbl>
## 1 Particle              4.816116e-07
## 2     Free              3.866798e-08

# Make a data frame to draw significance line in boxplot (visually calculated)
dat3 <- data.frame(a = c(1.1,1.1,1.9,1.9), b = c(-5.05,-5,-5,-5.05)) # WholePart vs WholeFree


poster_c <- ggplot(filter(metadata, norep_filter_name != "MOTEJ515"), 
       aes(y = log10(fracprod_per_cell), x = fraction)) +
  geom_jitter(size = 3, aes(fill = fraction, shape = season), width = 0.2) + 
  geom_boxplot(alpha = 0.5, outlier.shape = NA, aes( fill = fraction)) +
  scale_fill_manual(values = fraction_colors, guide = FALSE) +
  scale_shape_manual(values = season_shapes) +
  ylim(c(-8.5, -4.9)) + 
  ylab("Log10(Per-Capita Production) \n(μgC/cell/day)") +
  ##### Particle vs free per-cell production 
  geom_path(data = dat3, aes(x = a, y = b), linetype = 1, color = "gray40") +
  annotate("text", x=1.5, y=-4.95, size = 8, color = "gray40", label= paste("***")) +
  annotate("text", x=1.5, y=-5.15,  size = 3.5, color = "gray40",
           label= paste("p =", round(percellprod_wilcox$p.value, digits = 5))) +
  theme(legend.position = "none",
        axis.title.x = element_blank())


######## FIGURE 1
# legend
legend <- get_legend(poster_a)

row1_plots <- plot_grid(poster_a + theme(legend.position = "none"), poster_b, poster_c,
          labels = c("A", "B", "C"),
          ncol = 3, nrow = 1)

#fig_1 <- 
  plot_grid(row1_plots, season_legend,
                   ncol = 1, nrow = 2, 
                   rel_heights = c(1, 0.05))

Linear Regressions with Fraction Diversity

#################################### Bulk Measure Production #################################### 
################### Richness ################### 
summary(lm(frac_bacprod ~ mean, data = richness))  # All samples together

## 
## Call:
## lm(formula = frac_bacprod ~ mean, data = richness)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -15.057 -12.017  -5.654   9.256  46.605 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)  
## (Intercept) 18.501168   9.039579   2.047   0.0528 .
## mean        -0.003335   0.018911  -0.176   0.8616  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 15.54 on 22 degrees of freedom
## Multiple R-squared:  0.001412,   Adjusted R-squared:  -0.04398 
## F-statistic: 0.0311 on 1 and 22 DF,  p-value: 0.8616

# Particle-Associated 
lm_prod_vs_rich_PA <- lm(frac_bacprod ~ mean, data = filter(richness, fraction == "Particle"))
summary(lm_prod_vs_rich_PA)

## 
## Call:
## lm(formula = frac_bacprod ~ mean, data = filter(richness, fraction == 
##     "Particle"))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -8.1271 -2.8865 -0.4649  3.7537  9.8401 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)   
## (Intercept) -11.269531   5.717450  -1.971  0.07701 . 
## mean          0.038125   0.009868   3.863  0.00314 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.476 on 10 degrees of freedom
## Multiple R-squared:  0.5988, Adjusted R-squared:  0.5587 
## F-statistic: 14.93 on 1 and 10 DF,  p-value: 0.003143

# Cross validate particle-associated for a better estimate of the adjusted R-squared
cv_lm_prod_vs_rich_PA <- train(
      frac_bacprod ~ mean, data = filter(richness, fraction == "Particle"), 
      method ='lm', 
      trControl = trainControl(method ="repeatedcv", number = 3, repeats = 100), 
      tuneGrid = expand.grid(intercept = TRUE))

cv_lm_prod_vs_rich_PA$results   # Particle-Associated CV results 

##   intercept     RMSE  Rsquared      MAE   RMSESD RsquaredSD    MAESD
## 1      TRUE 6.714924 0.6410184 5.362125 2.214199  0.2966807 1.819534

summary(lm(frac_bacprod ~ mean, data = filter(richness, fraction == "Free"))) # Free Living 

## 
## Call:
## lm(formula = frac_bacprod ~ mean, data = filter(richness, fraction == 
##     "Free"))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -22.094 -11.112  -2.755   8.611  33.308 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  5.40793   21.62006   0.250    0.808
## mean         0.05511    0.06207   0.888    0.395
## 
## Residual standard error: 17.7 on 10 degrees of freedom
## Multiple R-squared:  0.07306,    Adjusted R-squared:  -0.01963 
## F-statistic: 0.7882 on 1 and 10 DF,  p-value: 0.3955

################### Shannon Entropy ################### 
summary(lm(frac_bacprod ~ mean, data = shannon))  # All samples together

## 
## Call:
## lm(formula = frac_bacprod ~ mean, data = shannon)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -14.659 -11.878  -5.666   9.231  46.593 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept) 16.63373   26.71652   0.623    0.540
## mean         0.08797    6.22969   0.014    0.989
## 
## Residual standard error: 15.55 on 22 degrees of freedom
## Multiple R-squared:  9.064e-06,  Adjusted R-squared:  -0.04545 
## F-statistic: 0.0001994 on 1 and 22 DF,  p-value: 0.9889

# Particle-Associated 
lm_prod_vs_shannon_PA <- lm(frac_bacprod ~ mean, data = filter(shannon, fraction == "Particle"))
summary(lm_prod_vs_shannon_PA)

## 
## Call:
## lm(formula = frac_bacprod ~ mean, data = filter(shannon, fraction == 
##     "Particle"))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.9022 -2.9150 -0.5875  1.6713 12.0544 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  -40.320     14.007  -2.878  0.01643 * 
## mean          11.089      3.068   3.614  0.00473 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.693 on 10 degrees of freedom
## Multiple R-squared:  0.5664, Adjusted R-squared:  0.5231 
## F-statistic: 13.06 on 1 and 10 DF,  p-value: 0.004734

# Cross validate particle-associated for a better estimate of the adjusted R-squared
cv_lm_prod_vs_shannon_PA <- train(
      frac_bacprod ~ mean, data = filter(shannon, fraction == "Particle"), 
      method ='lm', 
      trControl = trainControl(method ="repeatedcv", number = 3, repeats = 100), 
      tuneGrid = expand.grid(intercept = TRUE))

cv_lm_prod_vs_shannon_PA$results   # Particle-Associated CV results 

##   intercept     RMSE  Rsquared      MAE   RMSESD RsquaredSD    MAESD
## 1      TRUE 6.483534 0.6694728 5.092528 2.388309  0.2923481 1.952157

summary(lm(frac_bacprod ~ mean, data = filter(shannon, fraction == "Free"))) # Free Living 

## 
## Call:
## lm(formula = frac_bacprod ~ mean, data = filter(shannon, fraction == 
##     "Free"))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -22.004 -10.708  -3.738   6.632  37.129 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)   -13.37      73.87  -0.181    0.860
## mean            9.40      18.50   0.508    0.622
## 
## Residual standard error: 18.15 on 10 degrees of freedom
## Multiple R-squared:  0.02516,    Adjusted R-squared:  -0.07233 
## F-statistic: 0.2581 on 1 and 10 DF,  p-value: 0.6225

################### Inverse Simpson ################### 
summary(lm(frac_bacprod ~ mean, data = invsimps))  # All samples together

## 
## Call:
## lm(formula = frac_bacprod ~ mean, data = invsimps)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -15.269 -10.298  -4.916   5.866  46.452 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  12.1577     6.0225   2.019   0.0559 .
## mean          0.1629     0.1731   0.941   0.3570  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 15.25 on 22 degrees of freedom
## Multiple R-squared:  0.03868,    Adjusted R-squared:  -0.005019 
## F-statistic: 0.8851 on 1 and 22 DF,  p-value: 0.357

# Particle-Associated samples 
lm_prod_vs_invsimps_PA <- lm(frac_bacprod ~ mean, data = filter(invsimps, fraction == "Particle"))
summary(lm_prod_vs_invsimps_PA)

## 
## Call:
## lm(formula = frac_bacprod ~ mean, data = filter(invsimps, fraction == 
##     "Particle"))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.5589 -2.1093 -0.1969  0.8752  7.8282 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.45199    2.43866  -0.185  0.85666    
## mean         0.29344    0.05781   5.076  0.00048 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.571 on 10 degrees of freedom
## Multiple R-squared:  0.7204, Adjusted R-squared:  0.6925 
## F-statistic: 25.77 on 1 and 10 DF,  p-value: 0.0004804

# Cross validate particle-associated for a better estimate of the adjusted R-squared
cv_lm_prod_vs_invsimps_PA <- train(
      frac_bacprod ~ mean, data = filter(invsimps, fraction == "Particle"), 
      method ='lm', 
      trControl = trainControl(method ="repeatedcv", number = 3, repeats = 100), 
      tuneGrid = expand.grid(intercept = TRUE))

cv_lm_prod_vs_invsimps_PA$results       # Cross Validated PA results

##   intercept     RMSE  Rsquared      MAE   RMSESD RsquaredSD    MAESD
## 1      TRUE 5.181725 0.7603274 3.923763 2.163189  0.2516618 1.792041

#Free Living Samples 
summary(lm(frac_bacprod ~ mean, data = filter(invsimps, fraction == "Free")))

## 
## Call:
## lm(formula = frac_bacprod ~ mean, data = filter(invsimps, fraction == 
##     "Free"))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -20.765  -9.356  -4.445   6.116  36.017 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  11.0985    16.7052   0.664    0.521
## mean          0.5379     0.6598   0.815    0.434
## 
## Residual standard error: 17.8 on 10 degrees of freedom
## Multiple R-squared:  0.06233,    Adjusted R-squared:  -0.03143 
## F-statistic: 0.6648 on 1 and 10 DF,  p-value: 0.4339

################### Simpson's Evenness ################### 
summary(lm(frac_bacprod ~ mean, data = simpseven))  # All samples together

## 
## Call:
## lm(formula = frac_bacprod ~ mean, data = simpseven)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -17.000  -7.163  -3.815   3.392  45.845 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  -0.8524     9.2286  -0.092   0.9272  
## mean        274.1606   134.4253   2.040   0.0536 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 14.26 on 22 degrees of freedom
## Multiple R-squared:  0.159,  Adjusted R-squared:  0.1208 
## F-statistic:  4.16 on 1 and 22 DF,  p-value: 0.05359

# Particle-Associated 
lm_prod_vs_simpseven_PA <- lm(frac_bacprod ~ mean, data = filter(simpseven, fraction == "Particle"))
summary(lm_prod_vs_simpseven_PA)

## 
## Call:
## lm(formula = frac_bacprod ~ mean, data = filter(simpseven, fraction == 
##     "Particle"))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -10.973  -2.229  -1.086   1.356  12.380 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)   -3.829      4.539  -0.844  0.41865   
## mean         234.057     71.238   3.286  0.00821 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.995 on 10 degrees of freedom
## Multiple R-squared:  0.5191, Adjusted R-squared:  0.471 
## F-statistic: 10.79 on 1 and 10 DF,  p-value: 0.008211

# Cross validate particle-associated for a better estimate of the adjusted R-squared
cv_lm_prod_vs_simpseven_PA <- train(
      frac_bacprod ~ mean, data = filter(simpseven, fraction == "Particle"), 
      method ='lm', 
      trControl = trainControl(method ="repeatedcv", number = 3, repeats = 100), 
      tuneGrid = expand.grid(intercept = TRUE))

cv_lm_prod_vs_simpseven_PA$results   # Particle-Associated CV results 

##   intercept     RMSE  Rsquared      MAE   RMSESD RsquaredSD    MAESD
## 1      TRUE 6.190526 0.6495025 4.726163 2.854309   0.304804 2.126732

summary(lm(frac_bacprod ~ mean, data = filter(simpseven, fraction == "Free"))) # Free Living 

## 
## Call:
## lm(formula = frac_bacprod ~ mean, data = filter(simpseven, fraction == 
##     "Free"))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -19.389 -12.029  -3.306   4.668  39.946 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)    15.85      23.52   0.674    0.516
## mean          114.96     321.02   0.358    0.728
## 
## Residual standard error: 18.27 on 10 degrees of freedom
## Multiple R-squared:  0.01266,    Adjusted R-squared:  -0.08607 
## F-statistic: 0.1282 on 1 and 10 DF,  p-value: 0.7277

#################################### Per-Cell Production #################################### 
################### Richness ################### 
summary(lm(log10(fracprod_per_cell_noinf) ~ mean, data = richness))  # All samples together 

## 
## Call:
## lm(formula = log10(fracprod_per_cell_noinf) ~ mean, data = richness)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.8668 -0.2124  0.1045  0.2133  0.6473 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -8.4591277  0.2212633 -38.231  < 2e-16 ***
## mean         0.0028988  0.0004641   6.246  3.4e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3804 on 21 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.6501, Adjusted R-squared:  0.6334 
## F-statistic: 39.02 on 1 and 21 DF,  p-value: 3.395e-06

# Particle-Associated 
lm_percell_prod_vs_rich_PA <- lm(log10(fracprod_per_cell_noinf) ~ mean, data = filter(richness, fraction == "Particle"))
summary(lm_percell_prod_vs_rich_PA)

## 
## Call:
## lm(formula = log10(fracprod_per_cell_noinf) ~ mean, data = filter(richness, 
##     fraction == "Particle"))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.55347 -0.21545 -0.01066  0.12536  0.58830 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -8.0617648  0.3717671 -21.685 4.44e-09 ***
## mean         0.0023794  0.0006348   3.748  0.00457 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3506 on 9 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.6095, Adjusted R-squared:  0.5662 
## F-statistic: 14.05 on 1 and 9 DF,  p-value: 0.004567

# Cross validate particle-associated for a better estimate of the adjusted R-squared
cv_lm_percell_prod_vs_rich_PA <- train(
      log10(fracprod_per_cell_noinf) ~ mean, data = filter(richness, fraction == "Particle" & norep_filter_name != "MOTEJ515"), 
      method ='lm', 
      trControl = trainControl(method ="repeatedcv", number = 3, repeats = 100), 
      tuneGrid = expand.grid(intercept = TRUE))

cv_lm_percell_prod_vs_rich_PA$results      # Particle-Associated CV results 

##   intercept      RMSE  Rsquared       MAE   RMSESD RsquaredSD     MAESD
## 1      TRUE 0.4075388 0.6452733 0.3251794 0.142633  0.3137362 0.1268408

summary(lm(log10(fracprod_per_cell_noinf) ~ mean, data = filter(richness, fraction == "Free"))) # Free Living 

## 
## Call:
## lm(formula = log10(fracprod_per_cell_noinf) ~ mean, data = filter(richness, 
##     fraction == "Free"))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.71719 -0.13833  0.07155  0.23581  0.56949 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -8.135758   0.471253 -17.264 8.99e-09 ***
## mean         0.001657   0.001353   1.225    0.249    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3859 on 10 degrees of freedom
## Multiple R-squared:  0.1304, Adjusted R-squared:  0.04344 
## F-statistic: 1.499 on 1 and 10 DF,  p-value: 0.2488

################### Shannon Entropy ################### 
summary(lm(log10(fracprod_per_cell_noinf) ~ mean, data = shannon))  # All samples together 

## 
## Call:
## lm(formula = log10(fracprod_per_cell_noinf) ~ mean, data = shannon)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.03514 -0.15042 -0.03394  0.26568  0.82794 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -10.9036     0.7534 -14.472 2.14e-12 ***
## mean          0.8805     0.1763   4.993 6.09e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4348 on 21 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.5428, Adjusted R-squared:  0.521 
## F-statistic: 24.93 on 1 and 21 DF,  p-value: 6.09e-05

# Particle-Associated 
lm_percell_prod_vs_shannon_PA <- lm(log10(fracprod_per_cell_noinf) ~ mean, data = filter(shannon, fraction == "Particle"))
summary(lm_percell_prod_vs_shannon_PA)

## 
## Call:
## lm(formula = log10(fracprod_per_cell_noinf) ~ mean, data = filter(shannon, 
##     fraction == "Particle"))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.36686 -0.23571 -0.01330  0.03961  0.70312 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -9.8897     0.8820 -11.213 1.37e-06 ***
## mean          0.6993     0.1935   3.614  0.00562 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3584 on 9 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.5921, Adjusted R-squared:  0.5468 
## F-statistic: 13.06 on 1 and 9 DF,  p-value: 0.00562

# Cross validate particle-associated for a better estimate of the adjusted R-squared
cv_lm_percell_prod_vs_shannon_PA <- train(
      log10(fracprod_per_cell_noinf) ~ mean, data = filter(shannon, fraction == "Particle" & norep_filter_name != "MOTEJ515"), 
      method ='lm', 
      trControl = trainControl(method ="repeatedcv", number = 3, repeats = 100), 
      tuneGrid = expand.grid(intercept = TRUE))

cv_lm_percell_prod_vs_shannon_PA$results      # Particle-Associated CV results 

##   intercept     RMSE  Rsquared       MAE    RMSESD RsquaredSD     MAESD
## 1      TRUE 0.413426 0.6987126 0.3345857 0.1459863  0.3050087 0.1300826

summary(lm(log10(fracprod_per_cell_noinf) ~ mean, data = filter(shannon, fraction == "Free"))) # Free Living 

## 
## Call:
## lm(formula = log10(fracprod_per_cell_noinf) ~ mean, data = filter(shannon, 
##     fraction == "Free"))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.72444 -0.17785  0.08114  0.13946  0.67366 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -8.8666     1.6332  -5.429 0.000289 ***
## mean          0.3243     0.4091   0.793 0.446298    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4014 on 10 degrees of freedom
## Multiple R-squared:  0.05914,    Adjusted R-squared:  -0.03495 
## F-statistic: 0.6285 on 1 and 10 DF,  p-value: 0.4463

################### Inverse Simpson ################### 
summary(lm(log10(fracprod_per_cell_noinf) ~ mean, data = invsimps))  # All samples together 

## 
## Call:
## lm(formula = log10(fracprod_per_cell_noinf) ~ mean, data = invsimps)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.97302 -0.28199 -0.05285  0.32003  0.99088 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -7.86039    0.18153 -43.301  < 2e-16 ***
## mean         0.02355    0.00525   4.485 0.000204 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4596 on 21 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.4893, Adjusted R-squared:  0.4649 
## F-statistic: 20.12 on 1 and 21 DF,  p-value: 0.0002037

# Particle-Associated 
lm_percell_prod_vs_invsimps_PA <- lm(log10(fracprod_per_cell_noinf) ~ mean, data = filter(invsimps, fraction == "Particle"))
summary(lm_percell_prod_vs_invsimps_PA)

## 
## Call:
## lm(formula = log10(fracprod_per_cell_noinf) ~ mean, data = filter(invsimps, 
##     fraction == "Particle"))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.28651 -0.18384 -0.11125  0.07337  0.56444 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -7.360961   0.159490 -46.153 5.27e-12 ***
## mean         0.018087   0.003759   4.812 0.000958 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2969 on 9 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.7201, Adjusted R-squared:  0.689 
## F-statistic: 23.15 on 1 and 9 DF,  p-value: 0.0009581

# Cross validate particle-associated for a better estimate of the adjusted R-squared
cv_lm_percell_prod_vs_invsimps_PA <- train(
      log10(fracprod_per_cell_noinf) ~ mean, data = filter(invsimps, fraction == "Particle" & norep_filter_name != "MOTEJ515"), 
      method ='lm', 
      trControl = trainControl(method ="repeatedcv", number = 3, repeats = 100), 
      tuneGrid = expand.grid(intercept = TRUE))

cv_lm_percell_prod_vs_invsimps_PA$results      # Particle-Associated CV results 

##   intercept      RMSE  Rsquared       MAE    RMSESD RsquaredSD      MAESD
## 1      TRUE 0.3569485 0.7164722 0.2874183 0.1139091  0.3273592 0.09448645

summary(lm(log10(fracprod_per_cell_noinf) ~ mean, data = filter(invsimps, fraction == "Free"))) # Free Living 

## 
## Call:
## lm(formula = log10(fracprod_per_cell_noinf) ~ mean, data = filter(invsimps, 
##     fraction == "Free"))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.68269 -0.11512  0.01325  0.17742  0.63175 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -8.03513    0.35685 -22.517 6.72e-10 ***
## mean         0.01910    0.01409   1.355    0.205    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3803 on 10 degrees of freedom
## Multiple R-squared:  0.1551, Adjusted R-squared:  0.07064 
## F-statistic: 1.836 on 1 and 10 DF,  p-value: 0.2052

################### Simpson's Evenness ################### 
summary(lm(log10(fracprod_per_cell_noinf) ~ mean, data = simpseven))  # All samples together 

## 
## Call:
## lm(formula = log10(fracprod_per_cell_noinf) ~ mean, data = simpseven)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.06874 -0.41920  0.04553  0.34511  1.56565 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -7.4496     0.4116  -18.10 2.74e-14 ***
## mean          4.3470     6.0344    0.72    0.479    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6353 on 21 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.02411,    Adjusted R-squared:  -0.02236 
## F-statistic: 0.5189 on 1 and 21 DF,  p-value: 0.4792

# Particle-Associated 
lm_percell_prod_vs_simpseven_PA <- lm(log10(fracprod_per_cell_noinf) ~ mean, data = filter(simpseven, fraction == "Particle"))
summary(lm_percell_prod_vs_simpseven_PA)

## 
## Call:
## lm(formula = log10(fracprod_per_cell_noinf) ~ mean, data = filter(simpseven, 
##     fraction == "Particle"))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.4726 -0.2230 -0.1106  0.1248  0.6887 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -7.5941     0.2885 -26.324 7.96e-10 ***
## mean         15.1887     4.6341   3.278  0.00957 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3788 on 9 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.5441, Adjusted R-squared:  0.4935 
## F-statistic: 10.74 on 1 and 9 DF,  p-value: 0.009566

# Cross validate particle-associated for a better estimate of the adjusted R-squared
cv_lm_percell_prod_vs_simpseven_PA <- train(
      log10(fracprod_per_cell_noinf) ~ mean, data = filter(simpseven, fraction == "Particle" & norep_filter_name != "MOTEJ515"), 
      method ='lm', 
      trControl = trainControl(method ="repeatedcv", number = 3, repeats = 100), 
      tuneGrid = expand.grid(intercept = TRUE))

cv_lm_percell_prod_vs_simpseven_PA$results      # Particle-Associated CV results 

##   intercept      RMSE  Rsquared      MAE   RMSESD RsquaredSD     MAESD
## 1      TRUE 0.4501359 0.6702939 0.362118 0.144351   0.316542 0.1181178

summary(lm(log10(fracprod_per_cell_noinf) ~ mean, data = filter(simpseven, fraction == "Free"))) # Free Living 

## 
## Call:
## lm(formula = log10(fracprod_per_cell_noinf) ~ mean, data = filter(simpseven, 
##     fraction == "Free"))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.63112 -0.24216  0.01388  0.14672  0.77426 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -7.9274     0.5202 -15.240    3e-08 ***
## mean          4.9361     7.1008   0.695    0.503    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4041 on 10 degrees of freedom
## Multiple R-squared:  0.04609,    Adjusted R-squared:  -0.0493 
## F-statistic: 0.4832 on 1 and 10 DF,  p-value: 0.5028

R2 table

########## PUT A TABLE TOGETHER 
# Per Liter Production 
perliter_row1 <- c("Richness", "Per-Liter",
          round(summary(lm_prod_vs_rich_PA)$adj.r.squared, digits = 3),
          round(cv_lm_prod_vs_rich_PA$results$Rsquared, digits = 3), 
          round(cv_lm_prod_vs_rich_PA$results$RsquaredSD, digits = 3))

perliter_row2 <- c("Shannon_Entropy", "Per-Liter",
          round(summary(lm_prod_vs_shannon_PA)$adj.r.squared, digits = 3),
          round(cv_lm_prod_vs_shannon_PA$results$Rsquared, digits = 3), 
          round(cv_lm_prod_vs_shannon_PA$results$RsquaredSD, digits = 3))

perliter_row3 <- c("Inverse_Simpson", "Per-Liter",
          round(summary(lm_prod_vs_invsimps_PA)$adj.r.squared, digits = 3),
          round(cv_lm_prod_vs_invsimps_PA$results$Rsquared, digits = 3), 
          round(cv_lm_prod_vs_invsimps_PA$results$RsquaredSD, digits = 3))

perliter_row4 <- c("Simpsons_Evenness","Per-Liter",
          round(summary(lm_prod_vs_simpseven_PA)$adj.r.squared, digits = 3),
          round(cv_lm_prod_vs_simpseven_PA$results$Rsquared, digits = 3), 
          round(cv_lm_prod_vs_simpseven_PA$results$RsquaredSD, digits = 3))



# Per cell production 
percell_row1 <- c("Richness", "Per-Cell",
          round(summary(lm_percell_prod_vs_rich_PA)$adj.r.squared, digits = 3),
          round(cv_lm_percell_prod_vs_rich_PA$results$Rsquared, digits = 3), 
          round(cv_lm_percell_prod_vs_rich_PA$results$RsquaredSD, digits = 3))

percell_row2 <- c("Shannon_Entropy", "Per-Cell",
          round(summary(lm_percell_prod_vs_shannon_PA)$adj.r.squared, digits = 3),
          round(cv_lm_percell_prod_vs_shannon_PA$results$Rsquared, digits = 3), 
          round(cv_lm_percell_prod_vs_shannon_PA$results$RsquaredSD, digits = 3))

percell_row3 <- c("Inverse_Simpson", "Per-Cell",
          round(summary(lm_percell_prod_vs_invsimps_PA)$adj.r.squared, digits = 3),
          round(cv_lm_percell_prod_vs_invsimps_PA$results$Rsquared, digits = 3), 
          round(cv_lm_percell_prod_vs_invsimps_PA$results$RsquaredSD, digits = 3))

percell_row4 <- c("Simpsons_Evenness", "Per-Cell",
          round(summary(lm_percell_prod_vs_simpseven_PA)$adj.r.squared, digits = 3),
          round(cv_lm_percell_prod_vs_simpseven_PA$results$Rsquared, digits = 3), 
          round(cv_lm_percell_prod_vs_simpseven_PA$results$RsquaredSD, digits = 3))

r2_table <- as.data.frame(rbind(perliter_row1, perliter_row2, perliter_row3, perliter_row4,
                                percell_row1, percell_row2, percell_row3, percell_row4))
colnames(r2_table) <- c("Diversity_Metric", "Production_Measure","Adjusted_R2","CV_R2", "CV_R2_SD")
row.names(r2_table) = NULL

pander(r2_table,
               caption = "Supplemental Table 2: \n R2 estimates for heterotrophic production vs particle-associated diversity linear regressions.")

Supplemental Table 2: R2 estimates for heterotrophic production vs particle-associated diversity linear regressions.

Diversity_Metric	Production_Measure	Adjusted_R2	CV_R2	CV_R2_SD
Richness	Per-Liter	0.559	0.641	0.297
Shannon_Entropy	Per-Liter	0.523	0.669	0.292
Inverse_Simpson	Per-Liter	0.692	0.76	0.252
Simpsons_Evenness	Per-Liter	0.471	0.65	0.305
Richness	Per-Cell	0.566	0.645	0.314
Shannon_Entropy	Per-Cell	0.547	0.699	0.305
Inverse_Simpson	Per-Cell	0.689	0.716	0.327
Simpsons_Evenness	Per-Cell	0.493	0.67	0.317

Figure 2

######################################################### OBSERVED RICHNESS ######################################################### 
rich_fraction_wilcox <- wilcox.test(mean ~ fraction, data = richness)
rich_fraction_wilcox

## 
##  Wilcoxon rank sum test
## 
## data:  mean by fraction
## W = 129, p-value = 0.0004955
## alternative hypothesis: true location shift is not equal to 0

filter(richness) %>%
  group_by(fraction) %>%
  summarize(mean(mean), sd(mean))

## # A tibble: 2 x 3
##   fraction `mean(mean)` `sd(mean)`
##     <fctr>        <dbl>      <dbl>
## 1 Particle     556.7992  167.31404
## 2     Free     338.4242   85.98665

# Make a data frame to draw significance line in boxplot (visually calculated)
rich_nums <- data.frame(a = c(1.15,1.15,1.85,1.85), b = c(920, 930, 930, 920)) # WholePart vs WholeFree

rich_distribution_plot <- 
  ggplot(richness, aes(y = mean, x = fraction)) +
  scale_fill_manual(values = fraction_colors) +
  geom_jitter(size = 3, aes(fill = fraction, shape = season), width = 0.2) + 
  geom_boxplot(alpha = 0.5, outlier.shape = NA, aes(fill = fraction)) +
  scale_y_continuous(limits = c(150,980), breaks = seq(from = 0, to =925, by = 150)) + 
  xlab("Observed Richness") + xlab("Fraction") + 
  scale_shape_manual(values = season_shapes) +
  geom_path(data = rich_nums, aes(x = a, y = b), linetype = 1, color = "#424645") +
  annotate("text", x=1.5, y=980, size = 8, color = "#424645", label= paste("***"), angle = 90) +
  annotate("text", x=1.5, y=790, size = 4, color = "#424645",
           label= paste("p =", round(rich_fraction_wilcox$p.value, digits = 5))) +
  theme(legend.position = "none",# axis.title.y = element_blank(),
        axis.title.x = element_blank(), axis.text.x = element_blank()) +
  coord_flip()

################ Richness vs Community-wide (Per-Liter) Production 
## 1. Extract the R2 and p-value from the linear model: 
lm_lab_perliter_rich_PA <- paste("atop(R^2 ==", round(summary(lm_prod_vs_rich_PA)$adj.r.squared, digits = 2), ",",
             "p ==", round(unname(summary(lm_prod_vs_rich_PA)$coefficients[,4][2]), digits = 3), ")")

## 2. Plot Richness vs Community-wide (Per-Liter) Production 
prod_vs_rich_plot <-  
  ggplot(richness, aes(x=mean, y=frac_bacprod)) + 
  geom_errorbarh(aes(xmin = mean - sd, xmax = mean + sd, color = fraction), alpha = 0.7) + # X-axis errorbars
  geom_errorbar(aes(ymin = frac_bacprod - SD_frac_bacprod, ymax = frac_bacprod + SD_frac_bacprod, color = fraction)) +  # Y-axis errorbars
  geom_point(size = 3.5, color = "black", aes(fill = fraction, shape = season)) +
  scale_color_manual(values = fraction_colors) +
  scale_fill_manual(values = fraction_colors) +
  scale_shape_manual(values = season_shapes) +
  ylab("\n Community Production \n (μgC/L/day)") + 
  xlab("Observed Richness") +
  geom_smooth(data=filter(richness, fraction == "Particle"), method='lm', color = "#FF6600", fill = "#FF6600") + 
  scale_x_continuous(limits = c(150,980), breaks = seq(from = 0, to =925, by = 150)) + 
   # Add the R2 and p-value to the plot 
  annotate("text", x=790, y=45, label=lm_lab_perliter_rich_PA, parse = TRUE, color = "#FF6600", size = 4) +
  theme(legend.position = "none", axis.title.x = element_blank(), axis.text.x = element_blank()) 


################ Richness vs Per Capita (Per-Cell) Production 
## 1. Extract the R2 and p-value from the linear model: 
lm_lab_percell_rich_PA <- paste("atop(R^2 ==", round(summary(lm_percell_prod_vs_rich_PA)$adj.r.squared, digits = 2), ",",
             "p ==", round(unname(summary(lm_percell_prod_vs_rich_PA)$coefficients[,4][2]), digits = 3), ")")

## 2. Plot Richness vs Per Capita (Per-Cell) Production 
percell_prod_vs_rich_plot <- 
  ggplot(richness, aes(x=mean, y=log10(fracprod_per_cell_noinf))) + 
  geom_errorbarh(aes(xmin = mean - sd, xmax = mean + sd, color = fraction), alpha = 0.7) + # X-axis errorbars
  geom_point(size = 3.5,  color = "black", aes(fill = fraction, shape = season)) +
  scale_color_manual(values = fraction_colors) +
  scale_fill_manual(values = fraction_colors) +
  scale_shape_manual(values = season_shapes) +
  ylab("\n log10(Per-Capita Production)\n (μgC/cell/day)") +
  xlab("Observed Richness") +
  geom_smooth(data=filter(richness, fraction == "Particle"), method='lm', color = "#FF6600", fill = "#FF6600") + 
  scale_x_continuous(limits = c(150,980), breaks = seq(from = 0, to =925, by = 150)) + 
  #scale_y_continuous(limits = c(-8e-7,5e-6), breaks = seq(from = 0, to = 6e-7, by = 3e-7)) + 
  # Add the R2 and p-value to the plot 
  annotate("text", x=790, y=-7.5, label=lm_lab_percell_rich_PA, parse = TRUE, color = "#FF6600", size = 4) +
  theme(legend.title = element_blank(), legend.position ="bottom", 
        legend.text = element_text(size = 14))




######################################################### INVERSE SIMPSON ######################################################### 
invsimps_fraction_wilcox <- wilcox.test(mean ~ fraction, data = invsimps)
invsimps_fraction_wilcox

## 
##  Wilcoxon rank sum test
## 
## data:  mean by fraction
## W = 81, p-value = 0.6297
## alternative hypothesis: true location shift is not equal to 0

filter(invsimps) %>%
  group_by(fraction) %>%
  summarize(mean(mean), sd(mean))

## # A tibble: 2 x 3
##   fraction `mean(mean)` `sd(mean)`
##     <fctr>        <dbl>      <dbl>
## 1 Particle     35.47659  23.843325
## 2     Free     24.09219   8.136901

# Make a data frame to draw significance line in boxplot (visually calculated)
invsimps_nums <- data.frame(a = c(1.15,1.15,1.85,1.85), b = c(83,85,85,83)) # WholePart vs WholeFree

invsimps_distribution_plot <- 
  ggplot(invsimps, aes(y = mean, x = fraction)) +
  scale_fill_manual(values = fraction_colors) +
  scale_shape_manual(values = season_shapes) +
  geom_jitter(size = 3,  aes(fill = fraction, shape = season), width = 0.2) + 
  geom_boxplot(alpha = 0.5, outlier.shape = NA, aes(fill = fraction)) +
  scale_y_continuous(limits = c(0,85), breaks = seq(from = 0, to = 85, by = 20)) + 
  ylab("Inverse Simpson") + xlab("Fraction") + 
  geom_path(data = invsimps_nums, aes(x = a, y = b), linetype = 1, color = "#424645") +
  annotate("text", x=1.5, y=80, size = 4, color = "#424645", label= "NS") +
  theme(legend.position = "none", #axis.title.y = element_blank(),
        axis.title.x = element_blank(), axis.text.x = element_blank()) +
  coord_flip()


################ Inverse Simpson vs Community-wide (Per-Liter) Production 
## 1. Extract the R2 and p-value from the linear model: 
lm_lab_perliter_invsimps_PA <- paste("atop(R^2 ==", round(summary(lm_prod_vs_invsimps_PA)$adj.r.squared, digits = 2), ",",
             "p ==", round(unname(summary(lm_prod_vs_invsimps_PA)$coefficients[,4][2]), digits = 4), ")")

## 2. Plot Inverse Simpson vs Community-wide (Per-Liter) Production 
prod_vs_invsimps_plot <-  
  ggplot(invsimps, aes(x=mean, y=frac_bacprod)) + 
  geom_errorbarh(aes(xmin = mean - sd, xmax = mean + sd, color = fraction), alpha = 0.7) + # X-axis errorbars
  geom_errorbar(aes(ymin = frac_bacprod - SD_frac_bacprod, ymax = frac_bacprod + SD_frac_bacprod, color = fraction)) +  # Y-axis errorbars
  geom_point(size = 3.5, color = "black", aes(fill = fraction, shape = season)) +
  scale_color_manual(values = fraction_colors) +
  scale_fill_manual(values = fraction_colors) +
  scale_shape_manual(values = season_shapes) +
  ylab("Community Production \n (μgC/L/day)") + 
  xlab("Inverse Simpson") +
  geom_smooth(data=filter(invsimps, fraction == "Particle"), method='lm', color = "#FF6600", fill = "#FF6600") + 
  scale_x_continuous(limits = c(0,85), breaks = seq(from = 0, to = 85, by = 20)) + 
  # Add the R2 and p-value to the plot 
  annotate("text", x=70, y=45, label=lm_lab_perliter_invsimps_PA, parse = TRUE, color = "#FF6600", size = 4) +
  theme(legend.position = "none", axis.title.x = element_blank(), axis.text.x = element_blank()) 



################ Inverse Simpson vs Per Capitra (Per-Cell) Production 
## 1. Extract the R2 and p-value from the linear model: 
lm_lab_percell_invsimps_PA <- paste("atop(R^2 ==", round(summary(lm_percell_prod_vs_invsimps_PA)$adj.r.squared, digits = 2), ",",
             "p ==", round(unname(summary(lm_percell_prod_vs_invsimps_PA)$coefficients[,4][2]), digits = 4), ")")

## 2. Plot Inverse Simpson vs Per Capitra (Per-Cell) Production 
percell_prod_vs_invsimps_plot <- 
  ggplot(invsimps, aes(x=mean, y=log10(fracprod_per_cell_noinf), fill = fraction, color = fraction)) + 
  geom_errorbarh(aes(xmin = mean - sd, xmax = mean + sd, color = fraction), alpha = 0.7) + # X-axis errorbars
  geom_point(size = 3.5, color = "black", aes(shape = season)) +
  scale_color_manual(values = fraction_colors, guide = TRUE) +
  scale_fill_manual(values = fraction_colors, guide = FALSE) +
  scale_shape_manual(values = season_shapes) +
  ylab("log10(production/cell)\n (μgC/cell/day)") +
  xlab("Inverse Simpson") +
  geom_smooth(data=filter(invsimps, fraction == "Particle"), method='lm', color = "#FF6600", fill = "#FF6600") + 
  scale_x_continuous(limits = c(0,85), breaks = seq(from = 0, to = 85, by = 20)) + 
  #scale_y_continuous(limits = c(-8e-7,5e-6), breaks = seq(from = 0, to = 6e-7, by = 3e-7)) + 
  # Add the R2 and p-value to the plot 
  annotate("text", x=70, y=-7.5, label=lm_lab_percell_invsimps_PA, parse = TRUE, color = "#FF6600", size = 4) +
  guides(color = guide_legend(override.aes = list(shape = 15))) +
  theme(legend.title = element_blank(), legend.position ="bottom", 
        legend.text = element_text(size = 14))


# Plot Inverse simpson plots altogether
invsimps_plots <- plot_grid(invsimps_distribution_plot, prod_vs_invsimps_plot, percell_prod_vs_invsimps_plot,
          labels = c("B", "D", "F"), ncol = 1, nrow = 3,
          rel_heights = c(0.5, 0.8, 1.1),
          align = "v")


# All 3 richness plots together 
rich_plots <- plot_grid(rich_distribution_plot, prod_vs_rich_plot, 
                        percell_prod_vs_rich_plot + theme(legend.position = "none"),
          labels = c("A", "C", "E"), ncol = 1, nrow = 3,
          rel_heights = c(0.5, 0.8, 1),
          align = "v")


invsimps_plots_noyaxis <- 
  plot_grid(invsimps_distribution_plot + theme(axis.title.y = element_blank(), axis.text.y = element_blank()), 
            prod_vs_invsimps_plot + theme(axis.title.y = element_blank()), 
            percell_prod_vs_invsimps_plot + theme(axis.title.y = element_blank(), legend.position = "none"), 
          labels = c("B", "D", "F"), ncol = 1, nrow = 3,
          rel_heights = c(0.5, 0.8, 1),
          align = "v")

# Extract distribution plots of richness and inverse simpson
figure2_row1 <- plot_grid(rich_plots, invsimps_plots_noyaxis,
          ncol = 2, nrow = 1, rel_widths = c(1, 0.75),
          align = "h")


######## PLOT FIGURE 2
plot_grid(figure2_row1, season_legend,
                   ncol = 1, nrow = 2, 
                   rel_heights = c(1, 0.05))

Figure 3

Unweighted ses.mpd

######################################################### SES MPD DISTRIBUTION: UNWEIGHTED  



################ Figure 3A: Distribution of Particle and Free Unweighted Mean Pairwise distance
## 1. Is there a difference between particle and free Unweighted Mean Pairwise distance?
unweighted_fraction_wilcox <- wilcox.test(mpd.obs.z ~ fraction, data = unweighted_df)
unweighted_fraction_wilcox

## 
##  Wilcoxon rank sum test
## 
## data:  mpd.obs.z by fraction
## W = 30, p-value = 0.01449
## alternative hypothesis: true location shift is not equal to 0

# 2. What are the means of particle and free unweighted PD?
unweighted_df %>%
  group_by(fraction) %>%
  summarize(mean(mpd.obs.z))

## # A tibble: 2 x 2
##   fraction `mean(mpd.obs.z)`
##     <fctr>             <dbl>
## 1 Particle        -0.4971164
## 2     Free         0.4375532

## 3. Plot Unweighted Phylogenetic Diversity vs Richness  
# Make a data frame to draw significance line in boxplot (visually calculated)
dat4 <- data.frame(a = c(1.15,1.15,1.85,1.85), b = c(1.6,1.7,1.7,1.6)) # WholePart vs WholeFree

unweight_distribution_plot <- 
  ggplot(unweighted_df, aes(y = mpd.obs.z, x = fraction)) +
  scale_fill_manual(values = fraction_colors) +
  geom_jitter(size = 3, aes(fill = fraction, shape = season), width = 0.2) + 
  geom_boxplot(alpha = 0.5, outlier.shape = NA, aes(fill = fraction)) +
  scale_shape_manual(values = season_shapes) + 
  scale_y_continuous(limits = c(-1.5,1.75), breaks = seq(from = -1.5, to = 1.5, by = 0.5)) + 
  ylab("Standardized Effect Size \n Unweighted Mean Pairwise Dist") +
  xlab("Fraction") + 
  geom_hline(yintercept = 0, linetype = "dashed", size = 1.5) +
  ##### Particle vs free per-cell production 
  geom_path(data = dat4, aes(x = a, y = b), linetype = 1, color = "#424645") +
  annotate("text", x=1.5, y=1.35, size = 4, color = "#424645",
           label= paste("p =", round(unweighted_fraction_wilcox$p.value, digits = 2))) +
  theme(legend.position = "none",# axis.title.y = element_blank(),
        axis.title.x = element_blank(), axis.text.x = element_blank()) +
  coord_flip()


################ Figure 3B: Unweighted Phylogenetic Diversity vs Richness  
## 1. Is there a relationship between richness and Unweighted Mean Pairwise distance?
lm_richness_vs_unweight <- lm(mean ~ mpd.obs.z, data = unweighted_df)

  # Linear model results for particle-associated only 
  summary(lm(mean ~ mpd.obs.z, data = filter(unweighted_df, fraction == "Particle")))

## 
## Call:
## lm(formula = mean ~ mpd.obs.z, data = filter(unweighted_df, fraction == 
##     "Particle"))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -191.79 -112.16  -41.89   55.90  278.87 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   515.15      51.16  10.069 1.49e-06 ***
## mpd.obs.z     -83.78      49.91  -1.679    0.124    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 155 on 10 degrees of freedom
## Multiple R-squared:  0.2199, Adjusted R-squared:  0.1419 
## F-statistic: 2.818 on 1 and 10 DF,  p-value: 0.1241

  # Linear model results for free-living only 
  summary(lm(mean ~ mpd.obs.z, data = filter(unweighted_df, fraction == "Free")))

## 
## Call:
## lm(formula = mean ~ mpd.obs.z, data = filter(unweighted_df, fraction == 
##     "Free"))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -100.15  -66.22  -18.70   43.65  172.92 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  337.088     42.753   7.884 1.34e-05 ***
## mpd.obs.z      3.053     77.510   0.039    0.969    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 90.18 on 10 degrees of freedom
## Multiple R-squared:  0.0001551,  Adjusted R-squared:  -0.09983 
## F-statistic: 0.001551 on 1 and 10 DF,  p-value: 0.9694

## 2. Extract the R2 and p-value from the linear model: 
lm_lab_richness_vs_unweight<- paste("atop(R^2 ==", round(summary(lm_richness_vs_unweight)$adj.r.squared, digits = 2), ",",
             "p ==", round(unname(summary(lm_richness_vs_unweight)$coefficients[,4][2]), digits = 3), ")")

## 3. Plot Unweighted Phylogenetic Diversity vs Richness  
unweight_rich_vs_mpd_plot <- 
  ggplot(unweighted_df, aes(y = mean, x = mpd.obs.z)) +
  geom_errorbar(aes(ymin = mean - sd, ymax = mean + sd, color = fraction), alpha = 0.7) + # X-axis errorbars
  geom_vline(xintercept = 0, linetype = "dashed", size = 1.5) +
  geom_point(size = 3, aes(fill = fraction, shape = season)) + ylab("\n Observed  Richness") +
  scale_shape_manual(values = season_shapes) + 
  xlab("Standardized Effect Size \n Unweighted Mean Pairwise Dist") +
  scale_fill_manual(values = fraction_colors) +
  geom_smooth(method = "lm", color = "#424645", fill = "#424645", alpha = 0.3) + 
  scale_x_continuous(limits = c(-1.5,1.75), breaks = seq(from = -1.5, to = 1.5, by = 0.5)) + 
  # Add the R2 and p-value to the plot 
  annotate("text", x=0.75, y=750, label=lm_lab_richness_vs_unweight, parse = TRUE, color = "#424645", size = 4) +
  theme(legend.title = element_blank(), legend.position = "none", 
        axis.text.x = element_blank(), axis.title.x = element_blank())


################ Figure 3C: Unweighted Mean Pairwise distance vs Community-Wide Production
# Is there a relationship between Production and Unweighted Mean Pairwise distance?
summary(lm(frac_bacprod ~ mpd.obs.z, data = unweighted_df)) # NS 

## 
## Call:
## lm(formula = frac_bacprod ~ mpd.obs.z, data = unweighted_df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -20.416  -9.547  -3.120   8.193  44.030 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   17.125      3.102   5.520 1.51e-05 ***
## mpd.obs.z      3.901      3.768   1.035    0.312    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 15.19 on 22 degrees of freedom
## Multiple R-squared:  0.04645,    Adjusted R-squared:  0.003106 
## F-statistic: 1.072 on 1 and 22 DF,  p-value: 0.3118

  # Linear model results for particle-associated only 
  summary(lm(frac_bacprod ~ mpd.obs.z, data = filter(unweighted_df, fraction == "Particle")))

## 
## Call:
## lm(formula = frac_bacprod ~ mpd.obs.z, data = filter(unweighted_df, 
##     fraction == "Particle"))
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -9.145 -4.312 -1.824  3.402 19.527 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)    8.213      2.617   3.138   0.0105 *
## mpd.obs.z     -3.511      2.553  -1.376   0.1990  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.928 on 10 degrees of freedom
## Multiple R-squared:  0.1591, Adjusted R-squared:  0.07501 
## F-statistic: 1.892 on 1 and 10 DF,  p-value: 0.199

  # Linear model results for free-living only 
  summary(lm(frac_bacprod ~ mpd.obs.z, data = filter(unweighted_df, fraction == "Free")))

## 
## Call:
## lm(formula = frac_bacprod ~ mpd.obs.z, data = filter(unweighted_df, 
##     fraction == "Free"))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -21.996 -14.939   2.192   8.751  37.004 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)   18.191      8.398   2.166   0.0555 .
## mpd.obs.z     13.410     15.225   0.881   0.3991  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 17.71 on 10 degrees of freedom
## Multiple R-squared:  0.072,  Adjusted R-squared:  -0.0208 
## F-statistic: 0.7758 on 1 and 10 DF,  p-value: 0.3991

unweight_prod_vs_mpd_plot <- 
  ggplot(unweighted_df, aes(y = frac_bacprod, x = mpd.obs.z, fill = fraction)) +
  geom_errorbar(aes(ymin = frac_bacprod - SD_frac_bacprod, ymax = frac_bacprod + SD_frac_bacprod, color = fraction), alpha = 0.7) + # X-axis errorbars
  geom_vline(xintercept = 0, linetype = "dashed", size = 1.5) +
  geom_point(size = 3, aes(shape = season)) +
  scale_shape_manual(values = season_shapes) + 
  ylab("\n Community Production \n (μgC/L/day)") +
  xlab("Standardized Effect Size \n Unweighted Mean Pairwise Dist") +
  scale_fill_manual(values = fraction_colors) +
  scale_x_continuous(limits = c(-1.5,1.75), breaks = seq(from = -1.5, to = 1.5, by = 0.5)) + 
  theme(legend.title = element_blank(), legend.position = "none",
        axis.text.x = element_blank(), axis.title.x = element_blank())



################ Figure 3D: Unweighted Phylogenetic Diversity vs Per Capitra (Per-Cell) Production 
## 1. Run the linear model: Is there a relationship between PER-CELL PRODUCTION and Unweighted Mean Pairwise distance?
unweight_vs_percell <- lm(log10(fracprod_per_cell_noinf) ~ mpd.obs.z, data = unweighted_df)

  # Linear model results for particle-associated only 
  summary(lm(log10(fracprod_per_cell_noinf) ~ mpd.obs.z, data = filter(unweighted_df, fraction == "Particle")))

## 
## Call:
## lm(formula = log10(fracprod_per_cell_noinf) ~ mpd.obs.z, data = filter(unweighted_df, 
##     fraction == "Particle"))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.3795 -0.2419 -0.1651  0.0443  1.1815 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -6.9248     0.1711 -40.475 1.71e-11 ***
## mpd.obs.z    -0.3299     0.1627  -2.028   0.0732 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4649 on 9 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.3136, Adjusted R-squared:  0.2373 
## F-statistic: 4.112 on 1 and 9 DF,  p-value: 0.07319

  # Linear model results for free-living only 
  summary(lm(log10(fracprod_per_cell_noinf) ~ mpd.obs.z, data = filter(unweighted_df, fraction == "Free")))

## 
## Call:
## lm(formula = log10(fracprod_per_cell_noinf) ~ mpd.obs.z, data = filter(unweighted_df, 
##     fraction == "Free"))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.63978 -0.29782  0.07694  0.17610  0.76503 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -7.55622    0.19603 -38.545  3.3e-12 ***
## mpd.obs.z   -0.04303    0.35540  -0.121    0.906    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4135 on 10 degrees of freedom
## Multiple R-squared:  0.001463,   Adjusted R-squared:  -0.09839 
## F-statistic: 0.01466 on 1 and 10 DF,  p-value: 0.906

## 2. Extract the R2 and p-value from the linear model: 
lm_lab_percell_unweightPD <- paste("atop(R^2 ==", round(summary(unweight_vs_percell)$adj.r.squared, digits = 2), ",",
             "p ==", round(unname(summary(unweight_vs_percell)$coefficients[,4][2]), digits = 4), ")")

## 3. Plot Unweighted Phylogenetic Diversity vs Per Capitra (Per-Cell) Production 
unweight_percell_vs_mpd_plot <- 
  ggplot(unweighted_df, 
       aes(y = log10(fracprod_per_cell_noinf), x = mpd.obs.z)) +
  geom_vline(xintercept = 0, linetype = "dashed", size = 1.5) +
  geom_point(size = 3,  aes(fill = fraction, shape = season)) + 
  ylab("\n log10(Per-Capita Production)\n (μgC/cell/day)") +
  xlab("Unweighted Phylogenetic Diversity") +
  scale_fill_manual(values = fraction_colors) +
  scale_shape_manual(values = season_shapes) + 
  geom_smooth(method = "lm", color="#424645", fill = "#424645", alpha = 0.3) +
  scale_y_continuous(limits = c(-8.5,-5), breaks = c(-8, -7, -6)) + 
  scale_x_continuous(limits = c(-1.5,1.75), breaks = seq(from = -1.5, to = 1.5, by = 0.5)) + 
  # Add the R2 and p-value to the plot 
  annotate("text", x=0.75, y=-6, label=lm_lab_percell_unweightPD, parse = TRUE, color = "#424645", size = 4) +
  theme(legend.title = element_blank(), legend.position ="bottom", 
        legend.text = element_text(size = 14))





######## FIGURE 3
figure3_row1 <- plot_grid(unweight_distribution_plot, unweight_rich_vs_mpd_plot, unweight_prod_vs_mpd_plot, 
          unweight_percell_vs_mpd_plot + theme(legend.position = "none"), 
          labels = c("A", "B", "C", "D"), ncol = 1, nrow = 4,
          rel_heights = c(0.5, 1, 1, 1.2),
          align = "v")

plot_grid(figure3_row1, season_legend,
                   ncol = 1, nrow = 2, 
                   rel_heights = c(1, 0.05))

LASSO

Prepare the data for Ridge & Lasso Regressions

unweight_vals <- unweighted_df %>%
  dplyr::select(mpd.obs.z, norep_filter_name) %>%
  rename(Unweighted_PD = mpd.obs.z)
  
weighted_vals <- weighted_df %>%
  dplyr::select(mpd.obs.z, norep_filter_name) %>%
  rename(Weighted_PD = mpd.obs.z)


wide_all_divs <- otu_alphadiv %>%
  dplyr::select(norep_filter_name ,mean, measure) %>%
  spread(measure, mean)

lasso_data_df <- metadata_pca %>%
  left_join(wide_all_divs, by = "norep_filter_name") %>%
  left_join(unweight_vals, by = "norep_filter_name") %>%
  left_join(weighted_vals, by = "norep_filter_name") %>%
  dplyr::select(-c(project, year, Date, limnion, norep_water_name, dnaconcrep1, 
                   SD_perc_attached_bacprod, SE_total_bac_abund, SE_perc_attached_abund, SE_attached_bac))


### PARTICLE DATA 
lasso_data_df_particle <- lasso_data_df %>%
  filter(fraction == "Particle" ) %>%
  dplyr::select(-fraction)

lasso_data_df_particle_noprod <- lasso_data_df_particle %>%
  dplyr::select(-c(fracprod_per_cell, fracprod_per_cell_noinf, 
                   tot_bacprod, SD_tot_bacprod,SD_frac_bacprod, 
                   norep_filter_name, BGA_cellspermL, Turb_NTU,
                   lakesite, season, station)) 
  
percell_lasso_data_df_particle_noprod <- lasso_data_df_particle %>%
  dplyr::select(-c(fracprod_per_cell, frac_bacprod,
                   tot_bacprod, SD_tot_bacprod,SD_frac_bacprod, 
                   lakesite, season, norep_filter_name, station, BGA_cellspermL, Turb_NTU)) %>%
  dplyr::filter(Temp_C != 14.28) # Remove the missing row of data :( )


## ALL DATA 
all_dat_lasso_comm <- lasso_data_df %>%
    dplyr::select(-c(fracprod_per_cell, fracprod_per_cell_noinf, 
                   tot_bacprod, SD_tot_bacprod,SD_frac_bacprod, 
                   norep_filter_name, BGA_cellspermL, Turb_NTU,
                   lakesite, season, station, fraction)) 
  

all_dat_lasso_percapita <- lasso_data_df  %>%
  dplyr::select(-c(fracprod_per_cell, frac_bacprod,
                   tot_bacprod, SD_tot_bacprod,SD_frac_bacprod, fraction,
                   lakesite, season, norep_filter_name, station, BGA_cellspermL, Turb_NTU)) %>%
  dplyr::filter(Temp_C != 14.28) # Remove the missing row of data :( )

Perform Ridge and Lasso regression

Lasso: Bulk Community Production: PARTICLE

# Set the seed for reproducibility of the grid values 
set.seed(777)


################ PREPARE DATA ################ 
# Subset data needed and scale all o
scaled_comm_data <- 
  lasso_data_df_particle_noprod %>%    # Use data only for the particle samples  
  scale() %>%                          # Scale all of the variables to have mean =0 and sd = 1
  as.data.frame()                      # Make it a dataframe so that model.matrix function works (does not take a matrix)

# Set model parameters for community level data
## NOTE: there cannot be any data with NA
x = model.matrix(frac_bacprod ~ ., scaled_comm_data)[,-1]
y = scaled_comm_data$frac_bacprod
grid = 10^seq(10,-2,length = 100)

################ PREPARE TRAINING & TESTING DATA FOR CROSS VALIDATION ################ 
# Pull out test and training sets for cross validation
# We will use half the set to train the model and the 2nd half of the dataset to test the model 
train <- sample(1:nrow(x), nrow(x)/2)
test <- -train
y_test <- y[test]

################ RIDGE ################ 
# Run RIDGE regression with alpha = 0
ridge_divs_train <- glmnet(x[train,], y[train], alpha = 0, lambda = grid, thresh = 1e-12, standardize = TRUE)
par(mfrow = c(1,2)) 
plot(ridge_divs_train)

# Cross validation
cv_ridge_divs <- cv.glmnet(x[train,], y[train], alpha = 0)
plot(cv_ridge_divs)

best_ridge_lambda <- cv_ridge_divs$lambda.min
ridge_divs_pred <- predict(ridge_divs_train, s = best_ridge_lambda, newx = x[test,])
mean((ridge_divs_pred - y_test)^2) # Test MSE

## [1] 1.104789

## Run ridge on the entire dataset 
ridge_divs <- glmnet(x, y, alpha = 0, lambda = grid, standardize = TRUE)
par(mfrow = c(1,1))
plot(ridge_divs)

################ LASSO ################ 
# Run lasso regression with alpha = 1
lasso_divs_train <- glmnet(x[train,], y[train], alpha = 1, lambda = grid, standardize = TRUE)
par(mfrow = c(1,2))
plot(lasso_divs_train)

# Cross validation
cv_lasso_divs <- cv.glmnet(x[train,], y[train], alpha = 1)
plot(cv_lasso_divs)

best_lasso_lambda <- cv_lasso_divs$lambda.min
lasso_divs_pred <- predict(lasso_divs_train, s = best_lasso_lambda, newx = x[test,])
mean((lasso_divs_pred - y_test)^2)

## [1] 1.405385

## Run lasso on the entire dataset with the best lambda value 
lasso_divs <- glmnet(x, y, alpha = 1, lambda = grid, standardize = TRUE)
par(mfrow = c(1,1))
plot(lasso_divs)

plot(lasso_divs, xvar = "lambda", label = TRUE)

plot(lasso_divs, xvar = "dev", label = TRUE)

# What are the lasso coefficients? (Anything with a . is not selected by the model)
predict(lasso_divs, type = "coefficients", s = best_lasso_lambda)

## 35 x 1 sparse Matrix of class "dgCMatrix"
##                                   1
## (Intercept)           -1.499450e-17
## Sample_depth_m         .           
## Temp_C                 .           
## SpCond_uSpercm         .           
## TDS_mgperL             .           
## pH                     .           
## ORP_mV                 .           
## Chl_Lab_ugperL         .           
## Cl_mgperL              .           
## SO4_mgperL             .           
## NO3_mgperL             .           
## NH3_mgperL             .           
## TKN_mgperL             .           
## SRP_ugperL             .           
## TP_ugperL              .           
## Alk_mgperL             .           
## DO_mgperL              .           
## DO_percent             .           
## total_bac_abund        .           
## attached_bac           .           
## perc_attached_abund    .           
## perc_attached_bacprod  .           
## fraction_bac_abund     .           
## PC1                    .           
## PC2                    .           
## PC3                    .           
## PC4                    .           
## PC5                    .           
## PC6                    .           
## Richness               .           
## Shannon_Entropy        .           
## Inverse_Simpson        1.609355e-01
## Simpsons_Evenness      .           
## Unweighted_PD          .           
## Weighted_PD            .

The test MSE for ridge regression is 1.1047887 while the test MSE for lasso is 1.4053855.

Additionally, the lasso model uses Inverse Simpson as the best and only predictor of production!

Lasso: Per capita Production: PARTICLE

scaled_percapita_data <- percell_lasso_data_df_particle_noprod %>%
  dplyr::filter(!is.na(fracprod_per_cell_noinf)) %>%
  mutate(log10_percell = log10(fracprod_per_cell_noinf)) %>%
  dplyr::select(-c(fracprod_per_cell_noinf)) %>%
  scale() %>%
  as.data.frame()

set.seed(777)

# Set model parameters for community level data
## NOTE: there cannot be any data with NA
x = model.matrix(log10_percell ~ ., scaled_percapita_data)[,-1]
y = scaled_percapita_data$log10_percell
grid = 10^seq(10,-2,length = 100)

# Pull out test and training sets for cross validation
# We will use half the set to train the model and the 2nd half of the dataset to test the model 
train <- sample(1:nrow(x), nrow(x)/2)
test <- -train
y_test <- y[test]

################ RIDGE
# Run RIDGE regression with alpha = 0
ridge_divs_train <- glmnet(x[train,], y[train], alpha = 0, lambda = grid, thresh = 1e-12, standardize = FALSE)
par(mfrow = c(1,2)) 
plot(ridge_divs_train)

# Cross validation
cv_ridge_divs <- cv.glmnet(x[train,], y[train], alpha = 0)
plot(cv_ridge_divs)

best_ridge_lambda <- cv_ridge_divs$lambda.min
ridge_divs_pred <- predict(ridge_divs_train, s = best_ridge_lambda, newx = x[test,])
mean((ridge_divs_pred - y_test)^2) # Test MSE

## [1] 1.57874

## Run ridge on the entire dataset 
ridge_divs <- glmnet(x, y, alpha = 0, lambda = grid, standardize = FALSE)
par(mfrow = c(1,1))
plot(ridge_divs)

################ LASSO
# Run lasso regression with alpha = 1
lasso_divs_train <- glmnet(x[train,], y[train], alpha = 1, lambda = grid, standardize = FALSE)
par(mfrow = c(1,2))
plot(lasso_divs_train)

# Cross validation
cv_lasso_divs <- cv.glmnet(x[train,], y[train], alpha = 1)
plot(cv_lasso_divs)

best_lasso_lambda <- cv_lasso_divs$lambda.min
lasso_divs_pred <- predict(lasso_divs_train, s = best_lasso_lambda, newx = x[test,])
mean((lasso_divs_pred - y_test)^2)

## [1] 1.743068

## Run lasso on the entire dataset 
lasso_divs <- glmnet(x, y, alpha = 1, lambda = grid, standardize = FALSE)
par(mfrow = c(1,1))
plot(lasso_divs)

plot(lasso_divs, xvar = "lambda", label = TRUE)

plot(lasso_divs, xvar = "dev", label = TRUE)

# What are the lasso coefficients?
predict(lasso_divs, type = "coefficients", s = best_lasso_lambda)

## 35 x 1 sparse Matrix of class "dgCMatrix"
##                                   1
## (Intercept)           -7.117994e-16
## Sample_depth_m         .           
## Temp_C                -1.378715e-01
## SpCond_uSpercm         .           
## TDS_mgperL             .           
## pH                    -3.662284e-03
## ORP_mV                 .           
## Chl_Lab_ugperL         .           
## Cl_mgperL              .           
## SO4_mgperL             .           
## NO3_mgperL             .           
## NH3_mgperL             .           
## TKN_mgperL             .           
## SRP_ugperL             .           
## TP_ugperL              .           
## Alk_mgperL             .           
## DO_mgperL              .           
## DO_percent             .           
## total_bac_abund       -6.587419e-03
## attached_bac          -3.907083e-02
## perc_attached_abund    .           
## perc_attached_bacprod  .           
## fraction_bac_abund    -3.577572e-18
## PC1                    .           
## PC2                    .           
## PC3                    .           
## PC4                    .           
## PC5                    1.233827e-02
## PC6                    .           
## Richness               .           
## Shannon_Entropy        .           
## Inverse_Simpson        4.347245e-01
## Simpsons_Evenness      .           
## Unweighted_PD          .           
## Weighted_PD            .

Lasso: Community-wide Production: ALL SAMPLES

# Set the seed for reproducibility of the grid values 
set.seed(777)

################ PREPARE DATA ################ 
# Subset data needed and scale all o
scaled_comm_data_ALL <- 
  all_dat_lasso_percapita %>%    # Use data only for the particle samples  
  scale() %>%                          # Scale all of the variables to have mean =0 and sd = 1
  as.data.frame()                      # Make it a dataframe so that model.matrix function works (does not take a matrix)

# Set model parameters for community level data
## NOTE: there cannot be any data with NA
x = model.matrix(frac_bacprod ~ ., scaled_comm_data_ALL)[,-1]

## Error in eval(predvars, data, env): object 'frac_bacprod' not found

y = scaled_comm_data$frac_bacprod
grid = 10^seq(10,-2,length = 100)

################ PREPARE TRAINING & TESTING DATA FOR CROSS VALIDATION ################ 
# Pull out test and training sets for cross validation
# We will use half the set to train the model and the 2nd half of the dataset to test the model 
train <- sample(1:nrow(x), nrow(x)/2)
test <- -train
y_test <- y[test]


################ LASSO ################ 
# Run lasso regression with alpha = 1
lasso_divs_train <- glmnet(x[train,], y[train], alpha = 1, lambda = grid, standardize = TRUE)
par(mfrow = c(1,2))
plot(lasso_divs_train)

# Cross validation
cv_lasso_divs <- cv.glmnet(x[train,], y[train], alpha = 1)
plot(cv_lasso_divs)

best_lasso_lambda <- cv_lasso_divs$lambda.min
lasso_divs_pred <- predict(lasso_divs_train, s = best_lasso_lambda, newx = x[test,])
mean((lasso_divs_pred - y_test)^2)

## Error in mean((lasso_divs_pred - y_test)^2): dims [product 6] do not match the length of object [7]

## Run lasso on the entire dataset with the best lambda value 
lasso_divs <- glmnet(x, y, alpha = 1, lambda = grid, standardize = TRUE)

## Error in glmnet(x, y, alpha = 1, lambda = grid, standardize = TRUE): number of observations in y (12) not equal to the number of rows of x (11)

par(mfrow = c(1,1))
plot(lasso_divs)

plot(lasso_divs, xvar = "lambda", label = TRUE)

plot(lasso_divs, xvar = "dev", label = TRUE)

# What are the lasso coefficients? (Anything with a . is not selected by the model)
predict(lasso_divs, type = "coefficients", s = best_lasso_lambda)

## 35 x 1 sparse Matrix of class "dgCMatrix"
##                                   1
## (Intercept)           -5.352278e-16
## Sample_depth_m         .           
## Temp_C                -1.619342e-01
## SpCond_uSpercm         .           
## TDS_mgperL             .           
## pH                    -3.623885e-02
## ORP_mV                 .           
## Chl_Lab_ugperL         .           
## Cl_mgperL              .           
## SO4_mgperL             .           
## NO3_mgperL             .           
## NH3_mgperL             .           
## TKN_mgperL             .           
## SRP_ugperL             .           
## TP_ugperL              .           
## Alk_mgperL             .           
## DO_mgperL              .           
## DO_percent             .           
## total_bac_abund       -7.539204e-02
## attached_bac          -8.188182e-02
## perc_attached_abund    .           
## perc_attached_bacprod  .           
## fraction_bac_abund    -3.174675e-03
## PC1                    .           
## PC2                    .           
## PC3                    .           
## PC4                    1.594507e-02
## PC5                    1.035121e-01
## PC6                    .           
## Richness               .           
## Shannon_Entropy        .           
## Inverse_Simpson        4.266024e-01
## Simpsons_Evenness      .           
## Unweighted_PD          .           
## Weighted_PD           -1.690375e-02

Lasso: Per-capita Production: ALL SAMPLES

set.seed(777)

scaled_percapita_ALL <- all_dat_lasso_percapita %>%
  dplyr::filter(!is.na(fracprod_per_cell_noinf)) %>%
  mutate(log10_percell = log10(fracprod_per_cell_noinf)) %>%
  dplyr::select(-c(fracprod_per_cell_noinf, perc_attached_abund, fraction_bac_abund, attached_bac)) %>%
  scale() %>%
  as.data.frame()


# Set model parameters for community level data
## NOTE: there cannot be any data with NA
x = model.matrix(log10_percell ~ ., scaled_percapita_ALL)[,-1]
y = scaled_percapita_ALL$log10_percell
grid = 10^seq(10,-2,length = 100)

# Pull out test and training sets for cross validation
# We will use half the set to train the model and the 2nd half of the dataset to test the model 
train <- sample(1:nrow(x), nrow(x)/2)
test <- -train
y_test <- y[test]

################ RIDGE
# Run RIDGE regression with alpha = 0
ridge_divs_train <- glmnet(x[train,], y[train], alpha = 0, lambda = grid, thresh = 1e-12, standardize = TRUE)
par(mfrow = c(1,2)) 
plot(ridge_divs_train)

# Cross validation
cv_ridge_divs <- cv.glmnet(x[train,], y[train], alpha = 0)
plot(cv_ridge_divs)

best_ridge_lambda <- cv_ridge_divs$lambda.min
ridge_divs_pred <- predict(ridge_divs_train, s = best_ridge_lambda, newx = x[test,])
mean((ridge_divs_pred - y_test)^2) # Test MSE

## [1] 1.105553

## Run ridge on the entire dataset 
ridge_divs <- glmnet(x, y, alpha = 0, lambda = grid, standardize = TRUE)
par(mfrow = c(1,1))
plot(ridge_divs)

################ LASSO
# Run lasso regression with alpha = 1
lasso_divs_train <- glmnet(x[train,], y[train], alpha = 1, lambda = grid, standardize = TRUE)
par(mfrow = c(1,2))
plot(lasso_divs_train)

# Cross validation
cv_lasso_divs <- cv.glmnet(x[train,], y[train], alpha = 1)
plot(cv_lasso_divs)

best_lasso_lambda <- cv_lasso_divs$lambda.min
lasso_divs_pred <- predict(lasso_divs_train, s = best_lasso_lambda, newx = x[test,])
mean((lasso_divs_pred - y_test)^2)

## [1] 0.9141164

## Run lasso on the entire dataset 
lasso_divs <- glmnet(x, y, alpha = 1, lambda = grid, standardize = TRUE)
par(mfrow = c(1,1))
plot(lasso_divs)

plot(lasso_divs, xvar = "lambda", label = TRUE)

plot(lasso_divs, xvar = "dev", label = TRUE)

# What are the lasso coefficients?
predict(lasso_divs, type = "coefficients", s = best_lasso_lambda)

## 32 x 1 sparse Matrix of class "dgCMatrix"
##                                   1
## (Intercept)            5.315922e-16
## Sample_depth_m         .           
## Temp_C                 .           
## SpCond_uSpercm         .           
## TDS_mgperL             .           
## pH                    -9.932255e-02
## ORP_mV                 .           
## Chl_Lab_ugperL         .           
## Cl_mgperL              .           
## SO4_mgperL            -3.544802e-04
## NO3_mgperL             .           
## NH3_mgperL             .           
## TKN_mgperL             .           
## SRP_ugperL             .           
## TP_ugperL              .           
## Alk_mgperL             .           
## DO_mgperL              .           
## DO_percent             .           
## total_bac_abund        .           
## perc_attached_bacprod  .           
## PC1                    .           
## PC2                    .           
## PC3                    .           
## PC4                    .           
## PC5                    .           
## PC6                    .           
## Richness               4.317733e-01
## Shannon_Entropy        .           
## Inverse_Simpson        .           
## Simpsons_Evenness      .           
## Unweighted_PD         -9.098229e-02
## Weighted_PD            .

Supplemental Figures

Figure S1

Map of Muskegon Lake

# Load the ggmap package for plotting of Muskegon Lake
library(ggmap)

# Extract the coordinates for Muskegon Lake 
lakesite_coordinates <- read.csv("data/metadata/ML_GPS_Coordinates.csv") %>%
  dplyr::select(lakesite, Latitude, Longitude) %>%
  dplyr::rename(Lat = Latitude, Long = Longitude)

# Set the boundaries for the map 
map <- get_map(c(-86.35, 43.21, -86.235, 43.265),  #left/bottom/right/top
                 zoom = 13, maptype = "toner-background")

# Plot the Muskegon Lake Map 
ggmap(map) + xlab("Longitude") + ylab("Latitude") +
  theme(axis.text = element_text(size = 10, color = "black"),
        axis.title = element_text(size = 12, color = "black", face = "bold")) +
    geom_point(data = lakesite_coordinates, 
               aes(x=Long, y=Lat), color="red", size=4) +
    geom_text(data = lakesite_coordinates,
              aes(x=Long, y=Lat, label=lakesite), col="white", cex=4, vjust = 0.5,hjust = -0.25) 

Figure S2

Faith’s PD

surface_PAFL_otu_pruned_RAREFIED_rm2

## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 1891 taxa and 24 samples ]
## sample_data() Sample Data:       [ 24 samples by 70 sample variables ]
## tax_table()   Taxonomy Table:    [ 1891 taxa by 8 taxonomic ranks ]
## phy_tree()    Phylogenetic Tree: [ 1891 tips and 1889 internal nodes ]

## Calculate Faith's PD and species richness for sample 
faiths_pd_RAREFIED <- pd(comm_RAREFIED_rm2, phy_RAREFIED_rm2, include.root = FALSE)
faiths_pd_RAREFIED$norep_filter_name <- row.names(faiths_pd_RAREFIED)


# Join Faith's PD with the rest of the metadata 
meta_data_PD <- left_join(faiths_pd_RAREFIED, lasso_data_df, by = "norep_filter_name")

################ Figure S1: species richness and faith's PD
## 1. Run the linear model: Is there a correlation between species richness and faith's PD?
lm_PD_vs_SR <- lm(PD ~ SR, data = meta_data_PD)
summary(lm_PD_vs_SR)

## 
## Call:
## lm(formula = PD ~ SR, data = meta_data_PD)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.93991 -0.96092  0.03063  1.00013  2.60206 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 10.836131   0.986599   10.98 2.13e-10 ***
## SR           0.062600   0.002285   27.40  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.52 on 22 degrees of freedom
## Multiple R-squared:  0.9715, Adjusted R-squared:  0.9702 
## F-statistic: 750.9 on 1 and 22 DF,  p-value: < 2.2e-16

## 2. Extract the R2 and p-value from the linear model: 
lm_lab_PD_vs_SR<- paste("atop(R^2 ==", round(summary(lm_PD_vs_SR)$adj.r.squared, digits = 3), ",",
             "p ==", round(unname(summary(lm_PD_vs_SR)$coefficients[,4][2]), digits = 20), ")")

## 3. Plot Figure S1: species richness and faith's PD
ggplot(meta_data_PD, aes(y = PD, x = SR)) + 
  geom_point(size = 3) + 
  ylab("Faith's PD") + 
  xlab("Species Richness") +
  geom_smooth(method = "lm", color = "black") +
  # Add the R2 and p-value to the plot 
  annotate("text", x=600, y=35, label=lm_lab_PD_vs_SR, parse = TRUE, color = "black", size = 4) 

Figure S3

Correlations between Particle and Free of values in Figure 1

work_df <- metadata %>%
  dplyr::select(norep_filter_name, fraction, fraction_bac_abund, frac_bacprod, fracprod_per_cell_noinf) %>%
  mutate(norep_water_name = paste(substr(norep_filter_name, 1,4), substr(norep_filter_name, 6,9), sep = "")) %>%
  dplyr::select(-norep_filter_name)

part_work_df <- work_df %>%
  filter(fraction == "Particle") %>%
  rename(part_bacabund = fraction_bac_abund,
         part_prod = frac_bacprod, 
         part_percell_prod = fracprod_per_cell_noinf) %>%
  dplyr::select(-fraction)

free_work_df <- work_df %>%
  filter(fraction == "Free") %>%
  rename(free_bacabund = fraction_bac_abund,
         free_prod = frac_bacprod, 
         free_percell_prod = fracprod_per_cell_noinf) %>%
  dplyr::select(-fraction)

byfrac_df <- part_work_df %>%
  left_join(free_work_df, by = "norep_water_name")

byfrac_df$season <- substr(byfrac_df$norep_water_name, 5,5) # 7th letter = month sampled
byfrac_df$season <- as.character(byfrac_df$season)
byfrac_df$season <- ifelse(byfrac_df$season == "5", "Spring", 
                             ifelse(byfrac_df$season == "7", "Summer", 
                                    ifelse(byfrac_df$season == "9", "Fall",
                                           "NA")))
byfrac_df$season <- factor(byfrac_df$season, levels = c("Spring", "Summer", "Fall"))

summary(lm(log10(part_bacabund) ~ log10(free_bacabund), data = filter(byfrac_df, norep_water_name != "MOTE515")))

## 
## Call:
## lm(formula = log10(part_bacabund) ~ log10(free_bacabund), data = filter(byfrac_df, 
##     norep_water_name != "MOTE515"))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.52690 -0.08048  0.05903  0.19565  0.35255 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)
## (Intercept)            2.0900     3.2247   0.648    0.533
## log10(free_bacabund)   0.4264     0.5532   0.771    0.461
## 
## Residual standard error: 0.3111 on 9 degrees of freedom
## Multiple R-squared:  0.06194,    Adjusted R-squared:  -0.04229 
## F-statistic: 0.5943 on 1 and 9 DF,  p-value: 0.4605

plot1 <- ggplot(filter(byfrac_df, norep_water_name != "MOTE515"), 
       aes(x = log10(free_bacabund), y = log10(part_bacabund))) +
  xlab("Free") +  ylab("Particle") + 
  ggtitle("Log10(Bacterial Counts) \n (cells/mL)") +
  geom_point(size = 3, fill = "grey", aes(shape = season), width = 0.2) + 
  scale_shape_manual(values = season_shapes) +
  theme(legend.position = "bottom",
        legend.title = element_blank())




################ Community-Wide production correlation between particle and free
## 1. Run the linear regression 
lm_prod_corr <- lm(part_prod ~ free_prod, data = byfrac_df)

# 2. Extract the R2 and p-value from the linear model: 
lm_lab_perliter_PAvsFL <- paste("atop(R^2 ==", round(summary(lm_prod_corr)$adj.r.squared, digits = 3), ",",
            "p ==", round(unname(summary(lm_prod_corr)$coefficients[,4][2]), digits = 3), ")")

## 3. Plot Community production correlation between particle and free
plot2 <- ggplot(byfrac_df, aes(x = free_prod, y = part_prod)) +
  xlab("Free") +  ylab("Particle") + 
  ggtitle("Community Production \n(μgC/L/day)") +
  geom_point(size = 3, fill = "grey", aes(shape = season), width = 0.2) + 
  scale_shape_manual(values = season_shapes) +
  geom_smooth(method = "lm", color = "black")  + 
  # Add the R2 and p-value to the plot 
  annotate("text", x=20, y=30, label=lm_lab_perliter_PAvsFL, parse = TRUE, color = "black", size = 4) +
  theme(legend.position = "none")


################ Per-Capita production correlation between particle and free
## 1. Run the linear regression 
lm_percell_corr <- lm(log10(part_percell_prod) ~ log10(free_percell_prod), data = byfrac_df)

# 2. Extract the R2 and p-value from the linear model: 
lm_lab_percell_PAvsFL <- paste("atop(R^2 ==", round(summary(lm_percell_corr)$adj.r.squared, digits = 3), ",",
            "p ==", round(unname(summary(lm_percell_corr)$coefficients[,4][2]), digits = 3), ")")

## 3. Plot Per-Capita production correlation between particle and free
plot3 <- ggplot(byfrac_df,
       aes(x = log10(free_percell_prod), y = log10(part_percell_prod))) +  
  xlab("Free") +  ylab("Particle") + 
  ggtitle("log10(Per-Capita Production)\n (μgC/cell/day)") +
  geom_point(size = 3, fill = "grey", aes(shape = season), width = 0.2) + 
  scale_shape_manual(values = season_shapes) +
  geom_smooth(method = "lm", color = "black") + 
  # Add the R2 and p-value to the plot 
  annotate("text", y = -5.8, x=-7.9, label=lm_lab_percell_PAvsFL, parse = TRUE, color = "black", size = 4) +
  theme(legend.position = "none")


# Make the final multi-paneled Plot
# Extract the legend
legend_s1 <- get_legend(plot1)

# Make the 3 plots
top_row_S1 <- plot_grid(plot1 +theme(legend.position = "none"), plot2, plot3, 
          nrow = 1, ncol = 3, 
          labels = c("A", "B", "C"), 
          align = "h")

# Put the legend and 3 plots together
plot_grid(top_row_S1, legend_s1,
          rel_heights = c(1, 0.05),
          nrow = 2, ncol = 1)

Figure S4

Compare Diversity Estimates across Station and Season

long_div_df <- lasso_data_df %>%
  dplyr::select(norep_filter_name, lakesite:season, Richness:Weighted_PD) %>%
  gather(measure, mean, Richness:Weighted_PD) %>%
  mutate(measure = factor(measure, levels = 
                            c("Richness", "Shannon_Entropy", "Inverse_Simpson", "Simpsons_Evenness", "Unweighted_PD", "Weighted_PD")))


divs_PAFLA_plot_lakesite <- ggplot(long_div_df, aes(y = mean, x = lakesite, shape = lakesite)) +
  ylab("Mean Diversity Value\n By Lake Station") +
  facet_wrap(~measure, scale = "free_y", ncol = 6) +
    geom_point(size = 3, aes(fill = fraction), color = "black", position = position_jitterdodge()) +
  geom_boxplot(alpha = 0.5, outlier.shape = NA, color = "black", aes(fill = fraction)) +
  scale_fill_manual(values = fraction_colors) + 
  scale_shape_manual(values = lakesite_shapes) + 
  scale_color_manual(values = fraction_colors) + 
  theme(legend.position = "none", axis.title.x = element_blank(), 
        axis.text.x = element_text(angle = 30, vjust = 1, hjust = 1)) 

plot_A <- plot_grid(divs_PAFLA_plot_lakesite, lakesite_legend,
          nrow = 2, ncol = 1,  labels = c("A",""),
          rel_heights = c(1, 0.05))

## Error in plot_grid(divs_PAFLA_plot_lakesite, lakesite_legend, nrow = 2, : object 'lakesite_legend' not found

divs_PAFLA_plot_season <- ggplot(long_div_df, aes(y = mean, x = season, shape = season)) +
  ylab("Mean Diversity Value \n By Season") +
  facet_wrap(~measure, scale = "free_y", ncol = 6) +
   geom_point(size = 3, aes(fill = fraction), color = "black", position = position_jitterdodge()) +
  geom_boxplot(alpha = 0.5, outlier.shape = NA, color = "black", aes(fill = fraction)) +
  scale_fill_manual(values = fraction_colors) + 
  scale_shape_manual(values = season_shapes) + 
  scale_color_manual(values = fraction_colors) + 
  theme(legend.position = "none", axis.title.x = element_blank(), 
        axis.text.x = element_text(angle = 30, vjust = 1, hjust = 1)) 

plot_B <- plot_grid(divs_PAFLA_plot_season, season_legend,
          nrow = 2, ncol = 1, labels = c("B",""),
          rel_heights = c(1, 0.05))

###### PLOT FIGURE S6
plot_grid(plot_A, plot_B,
          nrow = 2, ncol = 1)

## Error in plot_grid(plot_A, plot_B, nrow = 2, ncol = 1): object 'plot_A' not found

Figure S5

Shannon and Simpson’s Evenness

######################################################### SHANNON ENTROPY ######################################################### 
shannon_fraction_wilcox <- wilcox.test(mean ~ fraction, data = shannon)
shannon_fraction_wilcox

## 
##  Wilcoxon rank sum test
## 
## data:  mean by fraction
## W = 119, p-value = 0.00556
## alternative hypothesis: true location shift is not equal to 0

filter(shannon) %>%
  group_by(fraction) %>%
  summarize(mean(mean), sd(mean))

## # A tibble: 2 x 3
##   fraction `mean(mean)` `sd(mean)`
##     <fctr>        <dbl>      <dbl>
## 1 Particle     4.534078  0.5594638
## 2     Free     3.982314  0.2958221

# Make a data frame to draw significance line in boxplot (visually calculated)
shannon_nums <- data.frame(a = c(1.15,1.15,1.85,1.85), b = c(5.8, 5.9, 5.9, 5.8)) # WholePart vs WholeFree

shannon_distribution_plot <- 
  ggplot(shannon, aes(y = mean, x = fraction)) +
  scale_fill_manual(values = fraction_colors) +
  geom_jitter(size = 3, aes(shape = season, fill = fraction), width = 0.2) + 
  geom_boxplot(alpha = 0.5, outlier.shape = NA, aes(fill = fraction)) +
  scale_y_continuous(limits = c(3.5,6.1), breaks = seq(from = 3, to =6, by = 0.5)) + 
  xlab("Shannon Entropy") +
  xlab("Fraction") + 
  scale_shape_manual(values = season_shapes) +
    # Add line and pval
  geom_path(data = shannon_nums, aes(x = a, y = b), linetype = 1, color = "#424645") +
  annotate("text", x=1.5, y=6.1, size = 8, color = "#424645", label= paste("**"), angle = 90) +
  annotate("text", x=1.5, y=5.5, size = 4, color = "#424645",
           label= paste("p =", round(shannon_fraction_wilcox$p.value, digits = 3))) +
  theme(legend.position = "none",# axis.title.y = element_blank(),
        axis.title.x = element_blank(), axis.text.x = element_blank()) +
  coord_flip()

################ Shannon vs Community-wide (Per-Liter) Production 
## 1. Extract the R2 and p-value from the linear model: 
lm_lab_perliter_shannon_PA <- paste("atop(R^2 ==", round(summary(lm_prod_vs_shannon_PA)$adj.r.squared, digits = 2), ",",
             "p ==", round(unname(summary(lm_prod_vs_shannon_PA)$coefficients[,4][2]), digits = 3), ")")

## 2. Plot Shannon vs Community-wide (Per-Liter) Production 
prod_vs_shannon_plot <-  
  ggplot(shannon, aes(x=mean, y=frac_bacprod)) + 
  geom_errorbarh(aes(xmin = mean - sd, xmax = mean + sd, color = fraction), alpha = 0.7) + # X-axis errorbars
  geom_errorbar(aes(ymin = frac_bacprod - SD_frac_bacprod, ymax = frac_bacprod + SD_frac_bacprod, color = fraction)) +  # Y-axis errorbars
  geom_point(size = 3.5, color = "black", aes(fill = fraction, shape = season)) +
  scale_color_manual(values = fraction_colors) +
  scale_shape_manual(values = season_shapes) +
  scale_fill_manual(values = fraction_colors) +
  ylab("\n Community Production \n (μgC/L/day)") + 
  scale_x_continuous(limits = c(3.5,6.1), breaks = seq(from = 3, to =6, by = 0.5)) + 
  xlab("Shannon Entropy") +
  geom_smooth(data=filter(shannon, fraction == "Particle"), method='lm', color = "#FF6600", fill = "#FF6600") + 
  # Add the R2 and p-value to the plot 
  annotate("text", x=5.25, y=45, label=lm_lab_perliter_shannon_PA, parse = TRUE, color = "#FF6600", size = 4) +
  theme(legend.position = "none", axis.title.x = element_blank(), axis.text.x = element_blank()) 


################ Shannon vs Per Capitra (Per-Cell) Production 
## 1. Extract the R2 and p-value from the linear model: 
lm_lab_percell_shannon_PA <- paste("atop(R^2 ==", round(summary(lm_percell_prod_vs_shannon_PA)$adj.r.squared, digits = 2), ",",
             "p ==", round(unname(summary(lm_percell_prod_vs_shannon_PA)$coefficients[,4][2]), digits = 3), ")")

## 2. Plot Shannon vs Per Capitra (Per-Cell) Production 
percell_prod_vs_shannon_plot <- 
  ggplot(shannon, aes(x=mean, y=log10(fracprod_per_cell_noinf))) + 
  geom_errorbarh(aes(xmin = mean - sd, xmax = mean + sd, color = fraction), alpha = 0.7) + # X-axis errorbars
  geom_point(size = 3.5, color = "black", aes(fill = fraction, shape = season)) +
  scale_color_manual(values = fraction_colors) +
  scale_fill_manual(values = fraction_colors) +
  scale_shape_manual(values = season_shapes) +
  ylab("\n log10(Per-Capita Production)\n (μgC/cell/day)") +
  xlab("Shannon Entropy") +
  geom_smooth(data=filter(shannon, fraction == "Particle"), method='lm', color = "#FF6600", fill = "#FF6600") + 
  scale_x_continuous(limits = c(3.5,6.1), breaks = seq(from = 3, to =6, by = 0.5)) + 
  #scale_y_continuous(limits = c(-8e-7,5e-6), breaks = seq(from = 0, to = 6e-7, by = 3e-7)) + 
  # Add the R2 and p-value to the plot 
  annotate("text", x=5.25, y=-7.5, label=lm_lab_percell_shannon_PA, parse = TRUE, color = "#FF6600", size = 4) +
  theme(legend.title = element_blank(), legend.position = "none")



shannon_plots <- plot_grid(shannon_distribution_plot, prod_vs_shannon_plot, percell_prod_vs_shannon_plot,
          labels = c("A", "C", "E"), ncol = 1, nrow = 3,
          rel_heights = c(0.5, 0.85, 1.1),
          align = "v")


######################################################### SIMPSON'S EVENNESS ######################################################### 
simpseven_fraction_wilcox <- wilcox.test(mean ~ fraction, data = simpseven)
simpseven_fraction_wilcox

## 
##  Wilcoxon rank sum test
## 
## data:  mean by fraction
## W = 55, p-value = 0.3474
## alternative hypothesis: true location shift is not equal to 0

filter(simpseven) %>%
  group_by(fraction) %>%
  summarize(mean(mean), sd(mean))

## # A tibble: 2 x 3
##   fraction `mean(mean)` `sd(mean)`
##     <fctr>        <dbl>      <dbl>
## 1 Particle   0.05890559 0.02537484
## 2     Free   0.07138806 0.01716014

# Make a data frame to draw significance line in boxplot (visually calculated)
simpseven_nums <- data.frame(a = c(1.15,1.15,1.85,1.85), b = c(0.14, 0.15, 0.15, 0.14)) # WholePart vs WholeFree

simpseven_distribution_plot <- 
  ggplot(simpseven, aes(y = mean, x = fraction)) +
  scale_fill_manual(values = fraction_colors) +
  geom_jitter(size = 3, aes(fill = fraction, shape = season), width = 0.2) + 
  geom_boxplot(alpha = 0.5, outlier.shape = NA, aes(fill = fraction)) +
  scale_y_continuous(limits = c(0,0.15), breaks = seq(from = 0, to = 0.15, by = 0.03)) + 
  ylab("Simpson's Evenness") +
  scale_shape_manual(values = season_shapes) +
  xlab("Fraction") + 
  geom_path(data = simpseven_nums, aes(x = a, y = b), linetype = 1, color = "#424645") +
  annotate("text", x=1.5, y=0.14, size = 4, color = "#424645", label= "NS") +
  theme(legend.position = "none",
        axis.title.x = element_blank(), axis.text.x = element_blank()) +
  coord_flip()




################ Simpson's Evenness vs Community-wide (Per-Liter) Production 
## 1. Extract the R2 and p-value from the linear model: 
lm_lab_perliter_simpseven_PA <- paste("atop(R^2 ==", round(summary(lm_prod_vs_simpseven_PA)$adj.r.squared, digits = 2), ",",
             "p ==", round(unname(summary(lm_prod_vs_simpseven_PA)$coefficients[,4][2]), digits = 3), ")")

## 2. Plot Simpson's Evenness vs Community-wide (Per-Liter) Production 
prod_vs_simpseven_plot <-  
  ggplot(simpseven, aes(x=mean, y=frac_bacprod)) + 
  geom_errorbarh(aes(xmin = mean - sd, xmax = mean + sd, color = fraction), alpha = 0.7) + # X-axis errorbars
  geom_errorbar(aes(ymin = frac_bacprod - SD_frac_bacprod, ymax = frac_bacprod + SD_frac_bacprod, color = fraction)) +  # Y-axis errorbars
  geom_point(size = 3.5, color = "black", aes(fill = fraction, shape = season)) +
  scale_color_manual(values = fraction_colors) +
  scale_shape_manual(values = season_shapes) +
  scale_fill_manual(values = fraction_colors) +
  ylab("\n Community Production \n (μgC/L/day)") + 
  ylab("Simpson's Evenness") +
  geom_smooth(data=filter(simpseven, fraction == "Particle"), method='lm', color = "#FF6600", fill = "#FF6600") + 
  scale_x_continuous(limits = c(0,0.15), breaks = seq(from = 0, to = 0.15, by = 0.03)) + 
  # Add the R2 and p-value to the plot 
  annotate("text", x=0.14, y=15, label=lm_lab_perliter_simpseven_PA, parse = TRUE, color = "#FF6600", size = 4) +
  theme(legend.position = "none", axis.title.x = element_blank(), axis.text.x = element_blank()) 



################ Simpson's Evenness vs Per Capitra (Per-Cell) Production 
## 1. Extract the R2 and p-value from the linear model: 
lm_lab_percell_simpseven_PA <- paste("atop(R^2 ==", round(summary(lm_percell_prod_vs_simpseven_PA)$adj.r.squared, digits = 2), ",",
             "p ==", round(unname(summary(lm_percell_prod_vs_simpseven_PA)$coefficients[,4][2]), digits = 3), ")")

## 2. Plot Simpson's Evenness vs Per Capitra (Per-Cell) Production 
percell_prod_vs_simpseven_plot <- 
  ggplot(simpseven, aes(x=mean, y=log10(fracprod_per_cell_noinf))) + 
  geom_errorbarh(aes(xmin = mean - sd, xmax = mean + sd, color = fraction), alpha = 0.7) + # X-axis errorbars
  geom_point(size = 3.5, color = "black", aes(fill = fraction, shape = season)) +
  scale_color_manual(values = fraction_colors) +
  scale_shape_manual(values = season_shapes) +
  scale_fill_manual(values = fraction_colors) +
  ylab("\n log10(Per-Capita Production)\n (μgC/cell/day)") +
  xlab("Simpson's Evenness") +
  geom_smooth(data=filter(simpseven, fraction == "Particle"), method='lm', color = "#FF6600", fill = "#FF6600") + 
  scale_x_continuous(limits = c(0,0.15), breaks = seq(from = 0, to = 0.15, by = 0.03)) + 
  # Add the R2 and p-value to the plot 
  annotate("text", x=0.14, y=-6.3, label=lm_lab_percell_simpseven_PA, parse = TRUE, color = "#FF6600", size = 4) +
  theme(legend.title = element_blank(), legend.position ="none")


simpseven_plots <- plot_grid(simpseven_distribution_plot, prod_vs_simpseven_plot, percell_prod_vs_simpseven_plot,
          labels = c("B", "D", "F"), ncol = 1, nrow = 3,
          rel_heights = c(0.5, 0.85, 1.1),
          align = "v")


#### Remove the axis labels for prettier multipaneled plotting
simpseven_plots_noyaxis <- 
  plot_grid(simpseven_distribution_plot + theme(axis.title.y = element_blank(), axis.text.y = element_blank()), 
            prod_vs_simpseven_plot + theme(axis.title.y = element_blank()), 
            percell_prod_vs_simpseven_plot + theme(axis.title.y = element_blank()), 
          labels = c("B", "D", "F"), ncol = 1, nrow = 3,
          rel_heights = c(0.5, 0.8, 1.1),
          align = "v")

figS2_row1 <- plot_grid(shannon_plots, simpseven_plots_noyaxis,
          ncol = 2, nrow = 1, rel_widths = c(1, 0.75),
          align = "h")

###### PLOT FIGURE S5
plot_grid(figS2_row1, season_legend,
          ncol = 1, nrow = 2, rel_heights = c(1, 0.05))

Figure S6

See analysis/OTU_Removal_Analysis* files on the GitHub Repo

Figure S7

Per capita production combined models

################ Figure S7A: Per-cell Heterotrophic Production vs All Observed Richness
## 1. Extract the R2 and p-value from the linear model: 
lm_lab_overall_rich_vs_percell <- 
  paste("atop(R^2 ==", round(summary(lm(log10(fracprod_per_cell_noinf) ~ mean, data = richness))$adj.r.squared, digits = 2), ",",
  "p ==", round(unname(summary(lm(log10(fracprod_per_cell_noinf) ~ mean, data = richness))$coefficients[,4][2]), digits = 6), ")")

## 2.  Figure S7A: Per-cell Heterotrophic Production vs All Observed Richness
plot_all_rich_percell <- 
  ggplot(richness, aes(x=mean, y=log10(fracprod_per_cell_noinf))) + 
  geom_errorbarh(aes(xmin = mean - sd, xmax = mean + sd), color = "grey", alpha = 0.7) + # X-axis errorbars
  geom_point(size = 3.5, aes(shape = season), color = "black", fill =  "grey") +
  scale_shape_manual(values = season_shapes) + 
  ylab("\n log10(Per-Capita Production)\n (μgC/cell/day)") + xlab("Observed Richness") +
  geom_smooth(method='lm', color = "#424645", fill = "#424645") + 
  scale_x_continuous(limits = c(150,950), breaks = seq(from = 0, to =925, by = 150)) + 
  scale_y_continuous(limits = c(-8.5, -5.4), breaks = seq(from = -9, to = -5, by = 1)) + 
  # Add the R2 and p-value to the plot 
  annotate("text", x=790, y=-8, label=lm_lab_overall_rich_vs_percell, parse = TRUE, color = "#424645", size = 4) +
  theme(legend.title = element_blank(), legend.position ="bottom", 
        legend.text = element_text(size = 14))


################ Figure S7B: Per-cell Heterotrophic Production vs All shannon
## 1. Extract the R2 and p-value from the linear model: 
lm_lab_overall_shannon_vs_percell <- 
  paste("atop(R^2 ==", round(summary(lm(log10(fracprod_per_cell_noinf) ~ mean, data = shannon))$adj.r.squared, digits = 2), ",",
  "p ==", round(unname(summary(lm(log10(fracprod_per_cell_noinf) ~ mean, data = shannon))$coefficients[,4][2]), digits = 5), ")")

## 2.  Figure S7B: Per-cell Heterotrophic Production vs All Shannon
plot_all_shannon_percell <- 
  ggplot(shannon, aes(x=mean, y=log10(fracprod_per_cell_noinf))) + 
  geom_errorbarh(aes(xmin = mean - sd, xmax = mean + sd), color = "grey", alpha = 0.7) + # X-axis errorbars
  geom_point(size = 3.5, aes(shape = season), color = "black", fill =  "grey") +
  scale_shape_manual(values = season_shapes) + 
  ylab("log10(production/cell)\n (μgC/cell/day)") + xlab("Shannon Entropy") +
  geom_smooth(method='lm', color = "#424645", fill = "#424645") + 
  scale_x_continuous(limits = c(3.5,6), breaks = seq(from = 3, to =6, by = 0.5)) + 
  scale_y_continuous(limits = c(-8.5, -5.4), breaks = seq(from = -9, to = -5, by = 1)) + 
  # Add the R2 and p-value to the plot 
  annotate("text", x=5.25, y=-8, label=lm_lab_overall_shannon_vs_percell, parse = TRUE, color = "#424645", size = 4) +
  theme(legend.title = element_blank(), legend.position ="bottom", 
        legend.text = element_text(size = 14))

################ Figure S7C: Per-cell Heterotrophic Production vs All Inverse Simpson
## 1. Extract the R2 and p-value from the linear model: 
lm_lab_overall_invsimps_vs_percell <- 
  paste("atop(R^2 ==", round(summary(lm(log10(fracprod_per_cell_noinf) ~ mean, data = invsimps))$adj.r.squared, digits = 2), ",",
  "p ==", round(unname(summary(lm(log10(fracprod_per_cell_noinf) ~ mean, data = invsimps))$coefficients[,4][2]), digits = 5), ")")

## 2.  Figure S7C: Per-cell Heterotrophic Production vs All Inverse Simpson
plot_all_invsimps_percell <- 
  ggplot(invsimps, aes(x=mean, y=log10(fracprod_per_cell_noinf))) + 
  geom_errorbarh(aes(xmin = mean - sd, xmax = mean + sd), color = "grey", alpha = 0.7) + # X-axis errorbars
  geom_point(size = 3.5, aes(shape = season), color = "black", fill =  "grey") +
  scale_shape_manual(values = season_shapes) + 
  ylab("log10(production/cell)\n (μgC/cell/day)") + xlab("Inverse Simpson") +
  geom_smooth(method='lm', color = "#424645", fill = "#424645") + 
  scale_x_continuous(limits = c(0,85), breaks = seq(from = 0, to = 85, by = 20)) + 
  scale_y_continuous(limits = c(-8.5, -5.4), breaks = seq(from = -9, to = -5, by = 1)) + 
  # Add the R2 and p-value to the plot 
  annotate("text", x=70, y=-8, label=lm_lab_overall_invsimps_vs_percell, parse = TRUE, color = "#424645", size = 4) +
  theme(legend.title = element_blank(), legend.position ="bottom", 
        legend.text = element_text(size = 14))

################ Figure S7D: Per-cell Heterotrophic Production vs All Simpson's Evenness
## 1. Extract the R2 and p-value from the linear model: 
lm_lab_overall_simpseven_vs_percell <- 
  paste("atop(R^2 ==", round(summary(lm(log10(fracprod_per_cell_noinf) ~ mean, data = simpseven))$adj.r.squared, digits = 2), ",",
  "p ==", round(unname(summary(lm(log10(fracprod_per_cell_noinf) ~ mean, data = simpseven))$coefficients[,4][2]), digits = 2), ")")

## 2.  Figure S7D: Per-cell Heterotrophic Production vs All Inverse Simpson
plot_all_simpseven_percell <- 
  ggplot(simpseven, aes(x=mean, y=log10(fracprod_per_cell_noinf))) + 
  geom_errorbarh(aes(xmin = mean - sd, xmax = mean + sd), color = "grey", alpha = 0.7) + # X-axis errorbars
  geom_point(size = 3.5, aes(shape = season), color = "black", fill =  "grey") +
  scale_shape_manual(values = season_shapes) + 
  ylab("log10(production/cell)\n (μgC/cell/day)") + xlab("Simpson's Evenness") +
  geom_smooth(method='lm', color = "#424645", fill = "#424645") + 
  scale_x_continuous(limits = c(0,0.15), breaks = seq(from = 0, to = 0.15, by = 0.03)) + 
  scale_y_continuous(limits = c(-8.5, -5.4), breaks = seq(from = -9, to = -5, by = 1)) + 
  # Add the R2 and p-value to the plot 
  annotate("text", x=0.12, y=-8, label=lm_lab_overall_simpseven_vs_percell, parse = TRUE, color = "#424645", size = 4) +
  theme(legend.title = element_blank(), legend.position ="bottom", 
        legend.text = element_text(size = 14))


figs3_row1 <- plot_grid(plot_all_rich_percell + theme(legend.position = "none"), 
          plot_all_shannon_percell + theme(axis.title.y = element_blank(), legend.position = "none"), 
          plot_all_invsimps_percell + theme(axis.title.y = element_blank(), legend.position = "none"), 
          plot_all_simpseven_percell + theme(axis.title.y = element_blank(), legend.position = "none"),
          align = "h", labels = c("A", "B", "C", "D"),
          rel_widths = c(1.05, 0.9, 0.9, 0.9),
          nrow = 1, ncol = 4)

######## FIGURE S7
# legend
legend_grey <- get_legend(plot_all_rich_percell)


plot_grid(figs3_row1, legend_grey,
                   ncol = 1, nrow = 2, 
                   rel_heights = c(1, 0.05))

Figure S8

Weighted ses.mpd

######################################################### SES MPD DISTRIBUTION: WEIGHTED 
weighted_fraction_wilcox <- wilcox.test(mpd.obs.z ~ fraction, data = weighted_df)
weighted_fraction_wilcox

## 
##  Wilcoxon rank sum test
## 
## data:  mpd.obs.z by fraction
## W = 81, p-value = 0.6297
## alternative hypothesis: true location shift is not equal to 0

filter(weighted_df) %>%
  group_by(fraction) %>%
  summarize(mean(mpd.obs.z))

## # A tibble: 2 x 2
##   fraction `mean(mpd.obs.z)`
##     <fctr>             <dbl>
## 1 Particle        -0.3607944
## 2     Free        -0.3854885

# Make a data frame to draw significance line in boxplot (visually calculated)
dat5 <- data.frame(a = c(1.15,1.15,1.85,1.85), b = c(1.6,1.7,1.7,1.6)) # WholePart vs WholeFree


weight_distribution_plot <- 
  ggplot(weighted_df, aes(y = mpd.obs.z, x = fraction)) +
  #scale_color_manual(values = fraction_colors) + 
  scale_fill_manual(values = fraction_colors) +
  geom_jitter(size = 3.5, aes(fill = fraction, shape = season), width = 0.2) + 
  geom_boxplot(alpha = 0.5, outlier.shape = NA, aes(fill = fraction)) +
  scale_shape_manual(values = season_shapes) + 
  scale_y_continuous(limits = c(-1.5,1.75), breaks = seq(from = -1.5, to = 1.5, by = 0.5)) + 
  ylab("Standardized Effect Size \n Unweighted Mean Pairwise Dist") +
  xlab("Fraction") + 
  geom_hline(yintercept = 0, linetype = "dashed", size = 1.5) +
  ##### Particle vs free per-cell production 
  geom_path(data = dat5, aes(x = a, y = b), linetype = 1, color = "#424645") +
  annotate("text", x=1.5, y=1.55, fontface = "bold",  size = 3.5, color = "#424645", label= "NS") +
  theme(legend.position = "none", #axis.title.y = element_blank(),
        axis.title.x = element_blank(), axis.text.x = element_blank()) +
  coord_flip()


# Is there a relationship between inverse simpson and Weighted Mean Pairwise distance?
#summary(lm(mean ~ mpd.obs.z, data = weighted_df)) # NS
  
summary(lm(mean ~ mpd.obs.z, data = filter(weighted_df, fraction == "Particle")))

## 
## Call:
## lm(formula = mean ~ mpd.obs.z, data = filter(weighted_df, fraction == 
##     "Particle"))
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -28.42 -19.15   0.40  12.95  40.67 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)    23.58      12.61   1.870   0.0911 .
## mpd.obs.z     -32.98      29.43  -1.121   0.2886  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 23.57 on 10 degrees of freedom
## Multiple R-squared:  0.1116, Adjusted R-squared:  0.02276 
## F-statistic: 1.256 on 1 and 10 DF,  p-value: 0.2886

summary(lm(mean ~ mpd.obs.z, data = filter(weighted_df, fraction == "Free")))

## 
## Call:
## lm(formula = mean ~ mpd.obs.z, data = filter(weighted_df, fraction == 
##     "Free"))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -12.0597  -4.3465  -0.8155   3.8832  18.4294 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   22.187      4.720   4.701  0.00084 ***
## mpd.obs.z     -4.942     10.486  -0.471  0.64753    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.441 on 10 degrees of freedom
## Multiple R-squared:  0.02173,    Adjusted R-squared:  -0.0761 
## F-statistic: 0.2221 on 1 and 10 DF,  p-value: 0.6475

weight_invsimps_vs_mpd_plot <- 
  ggplot(weighted_df, aes(y = mean, x = mpd.obs.z, fill = fraction)) +
  geom_errorbar(aes(ymin = mean - sd, ymax = mean + sd, color = fraction), alpha = 0.7) + # X-axis errorbars
  geom_vline(xintercept = 0, linetype = "dashed", size = 1.5) +
  geom_point(size = 3.5, aes(shape = season)) + ylab("Inverse Simpson") +
  scale_shape_manual(values = season_shapes) + 
  #xlab("Standardized Effect Size \n Weighted Mean Pairwise Dist") +
  scale_fill_manual(values = fraction_colors) +
  scale_x_continuous(limits = c(-1.5,1.75), breaks = seq(from = -1.5, to = 1.5, by = 0.5)) + 
  theme(legend.title = element_blank(), legend.position = "none",
        axis.text.x = element_blank(), axis.title.x = element_blank())




# Is there a relationship between PER-LITER PRODUCTION and WEIGHTED Mean Pairwise distance?
summary(lm(frac_bacprod ~ mpd.obs.z, data = weighted_df))

## 
## Call:
## lm(formula = frac_bacprod ~ mpd.obs.z, data = weighted_df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -18.030 -10.236  -2.842   5.937  38.409 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)    9.008      5.664   1.590    0.126
## mpd.obs.z    -21.439     12.889  -1.663    0.110
## 
## Residual standard error: 14.66 on 22 degrees of freedom
## Multiple R-squared:  0.1117, Adjusted R-squared:  0.07133 
## F-statistic: 2.767 on 1 and 22 DF,  p-value: 0.1104

summary(lm(frac_bacprod ~ mpd.obs.z, data = filter(weighted_df, fraction == "Particle")))

## 
## Call:
## lm(formula = frac_bacprod ~ mpd.obs.z, data = filter(weighted_df, 
##     fraction == "Particle"))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -10.240  -4.575  -1.717   4.503  15.352 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)    4.391      4.126   1.064    0.312
## mpd.obs.z    -15.432      9.627  -1.603    0.140
## 
## Residual standard error: 7.711 on 10 degrees of freedom
## Multiple R-squared:  0.2044, Adjusted R-squared:  0.1249 
## F-statistic: 2.569 on 1 and 10 DF,  p-value: 0.14

summary(lm(frac_bacprod ~ mpd.obs.z, data = filter(weighted_df, fraction == "Free")))

## 
## Call:
## lm(formula = frac_bacprod ~ mpd.obs.z, data = filter(weighted_df, 
##     fraction == "Free"))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -19.790 -12.275  -3.481   7.360  30.573 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)   14.697      9.683   1.518    0.160
## mpd.obs.z    -24.285     21.512  -1.129    0.285
## 
## Residual standard error: 17.32 on 10 degrees of freedom
## Multiple R-squared:  0.113,  Adjusted R-squared:  0.02434 
## F-statistic: 1.274 on 1 and 10 DF,  p-value: 0.2853

weight_prod_vs_mpd_plot <- 
  ggplot(weighted_df, aes(y = frac_bacprod, x = mpd.obs.z, fill = fraction)) +
  geom_errorbar(aes(ymin = frac_bacprod - SD_frac_bacprod, ymax = frac_bacprod + SD_frac_bacprod, color = fraction), alpha = 0.7) + # X-axis errorbars
  geom_vline(xintercept = 0, linetype = "dashed", size = 1.5) +
  geom_point(size = 3.5, aes(shape = season)) + 
  scale_shape_manual(values = season_shapes) + 
  ylab("Heterotrophic Production \n (μgC/L/day)") +
  xlab("Standardized Effect Size \n Weighted Mean Pairwise Dist") +
  scale_fill_manual(values = fraction_colors) +
  scale_x_continuous(limits = c(-1.5,1.75), breaks = seq(from = -1.5, to = 1.5, by = 0.5)) + 
  theme(legend.title = element_blank(), legend.position = "none",
        axis.text.x = element_blank(), axis.title.x = element_blank())






################ Weighted Phylogenetic Diversity vs Per Capita (Per-Cell) Production 
# Is there a relationship between PER-CELL PRODUCTION and WEIGHTED Mean Pairwise distance?
summary(lm(log10(fracprod_per_cell_noinf) ~ mpd.obs.z, data = weighted_df))

## 
## Call:
## lm(formula = log10(fracprod_per_cell_noinf) ~ mpd.obs.z, data = weighted_df)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.94765 -0.40616 -0.03619  0.31885  1.39101 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -7.5352     0.2442 -30.860   <2e-16 ***
## mpd.obs.z    -0.9519     0.5446  -1.748   0.0951 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6009 on 21 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.127,  Adjusted R-squared:  0.08544 
## F-statistic: 3.055 on 1 and 21 DF,  p-value: 0.09508

summary(lm(log10(fracprod_per_cell_noinf) ~ mpd.obs.z, data = filter(weighted_df, fraction == "Particle")))

## 
## Call:
## lm(formula = log10(fracprod_per_cell_noinf) ~ mpd.obs.z, data = filter(weighted_df, 
##     fraction == "Particle"))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.65329 -0.32632 -0.03486  0.22224  0.95307 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -7.0845     0.3015 -23.496 2.18e-09 ***
## mpd.obs.z    -0.9337     0.6753  -1.383      0.2    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5096 on 9 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.1752, Adjusted R-squared:  0.08357 
## F-statistic: 1.912 on 1 and 9 DF,  p-value: 0.2001

lm_percell_free_mpd <- lm(log10(fracprod_per_cell_noinf) ~ mpd.obs.z, data = filter(weighted_df, fraction == "Free"))
summary(lm_percell_free_mpd)

## 
## Call:
## lm(formula = log10(fracprod_per_cell_noinf) ~ mpd.obs.z, data = filter(weighted_df, 
##     fraction == "Free"))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.54189 -0.16159 -0.00204  0.20070  0.44618 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -7.9519     0.1848 -43.019 1.11e-12 ***
## mpd.obs.z    -0.9777     0.4107  -2.381   0.0386 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3306 on 10 degrees of freedom
## Multiple R-squared:  0.3617, Adjusted R-squared:  0.2979 
## F-statistic: 5.667 on 1 and 10 DF,  p-value: 0.03857

## 1. Extract the R2 and p-value from the linear model: 
lm_lab_percell_free_mpd <- paste("atop(R^2 ==", round(summary(lm_percell_free_mpd)$adj.r.squared, digits = 2), ",",
             "p ==", round(unname(summary(lm_percell_free_mpd)$coefficients[,4][2]), digits = 3), ")")

## 2. Plot Weighted Pylogenetic Diersity vs Per Capita (Per-Cell) Production 
weight_percell_vs_mpd_plot <- 
  ggplot(weighted_df, 
       aes(y = log10(fracprod_per_cell_noinf), x = mpd.obs.z, fill = fraction)) +
  geom_vline(xintercept = 0, linetype = "dashed", size = 1.5) +
  geom_point(size = 3.5, aes(shape = season)) + 
  scale_shape_manual(values = season_shapes) + 
  ylab("log10(Per-Capita Production)\n (μgC/cell/day)") +
  xlab("Weighted Phylogenetic Diversity") +
  scale_fill_manual(values = fraction_colors) +
  geom_smooth(method = "lm", data = filter(weighted_df, fraction == "Free"), color = "skyblue", fill = "skyblue", alpha = 0.3) +
  scale_x_continuous(limits = c(-1.5,1.75), breaks = seq(from = -1.5, to = 1.5, by = 0.5)) + 
  scale_y_continuous(limits = c(-8.5,-5), breaks = c(-8, -7, -6)) + 
  # Add the R2 and p-value to the plot 
  annotate("text", x=0.65, y=-7.9, label=lm_lab_percell_free_mpd, parse = TRUE, color = "deepskyblue3", size = 4) +
  theme(legend.title = element_blank(), legend.position ="none", 
        legend.text = element_text(size = 14))


row1_weighted_PD_plot <- plot_grid(weight_distribution_plot, weight_invsimps_vs_mpd_plot, weight_prod_vs_mpd_plot, weight_percell_vs_mpd_plot, 
          labels = c("A", "B", "C", "D"), ncol = 1, nrow = 4,
          rel_heights = c(0.5, 0.8, 0.8, 1.2),
          align = "v")

### PLOT FIGURE S8
plot_grid(row1_weighted_PD_plot, season_legend,
          ncol = 1, nrow = 2,
          rel_heights = c(1, 0.05))

Figure S9

Randomized Richness

#########################################################  Subset only richness data 
### These are the richness values for the fake samples 
#rich_stats <- filter(otu_alphadiv, measure == "Richness") %>%
#  dplyr::select(1:2) %>%
#  rename(mean_richness = mean) %>%
#  mutate(sample = paste("Sample_", seq(1:nrow(filter(otu_alphadiv, measure == "Richness"))), sep = ""),
#         mean_richness = matround(mean_richness))

## Pick OTUs to match these richness values 

  # List the otus from ALL samples 
#  all_otus <- taxa_names(surface_PAFL_otu_pruned_rm2)
  
  # Obtain the OTU table from the phyloseq object
#  otutab <- otu_table(surface_PAFL_otu_pruned_rm2)
  # Make all the counts to be 0
#  otutab_newvals <- apply(otutab, c(1, 2), function(x) 0)

  # Stop if things are wrong 
#  stopifnot(colnames(otutab_newvals) == all_otus)                       # Do the OTU names match?
#  stopifnot(rownames(otutab_newvals) == rich_stats$norep_filter_name)   # Do the sample names match?
  
# Make it reproducible!   
#set.seed(777)
  
#for (row in 1:nrow(rich_stats)) {
  
  # Pick the richness value 
#  rich_val <- rich_stats[row, 2]  
  
  # Sample the OTUs to represent that richness value 
#  col_index <- sample(ncol(otutab_newvals), rich_val, replace = FALSE, prob = NULL)
  
  # make all other columns 0
#  otutab_newvals[row, col_index] <- 1

#}


## Calculate the tree for those randomized samples 
# create a new phyloseq object 
#random_physeq_presab_raw <- phyloseq(otu_table(otutab_newvals, taxa_are_rows = FALSE), 
#                                 tax_table(surface_PAFL_otu_pruned_rm2), sample_data(surface_PAFL_otu_pruned_rm2))
#random_physeq_presab_raw

# Remove taxa that are 0!
#random_physeq_presab_pruned <- prune_taxa(taxa_sums(random_physeq_presab_raw) > 0, random_physeq_presab_raw) 
#random_physeq_presab_pruned

# Calculate tree 
# Obtain the OTU names that were retained in the dataset
#tax <- data.frame(tax_table(random_physeq_presab_pruned))
#vector_of_otus <- as.vector(tax$OTU)

# Write out the file for processing/fasttree
# write(vector_of_otus, file = "data/PhyloTree/randomized/random_physeq_presab_pruned_OTUnames.txt", ncolumns = 1, append = FALSE, sep = "\n")

# Read in the tree
randomized_tree <- read.tree(file = "data/PhyloTree/randomized/newick_tree_randomized_rmN.tre")
  
  #random_physeq_presab_pruned_tree <- merge_phyloseq(random_physeq_presab_pruned, randomized_tree)
  #random_physeq_presab_pruned_tree
#save(list="random_physeq_presab_pruned_tree", file=paste0("data/PhyloTree/randomized/random_physeq_presab_pruned_tree.RData")) 

load("data/PhyloTree/randomized/random_physeq_presab_pruned_tree.RData")
random_physeq_presab_pruned_tree

## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 2911 taxa and 24 samples ]
## sample_data() Sample Data:       [ 24 samples by 44 sample variables ]
## tax_table()   Taxonomy Table:    [ 2911 taxa by 8 taxonomic ranks ]
## phy_tree()    Phylogenetic Tree: [ 2911 tips and 2909 internal nodes ]

# First force the OTU 
randomized_otu <- matrix(otu_table(random_physeq_presab_pruned_tree), nrow = nrow(otu_table(random_physeq_presab_pruned_tree)))
rownames(randomized_otu) <- sample_names(random_physeq_presab_pruned_tree)
colnames(randomized_otu) <- taxa_names(random_physeq_presab_pruned_tree)
    
  
## Calculate input for SES_MPD  
# Convert the presence/absence data to standardized abundanced  vegan function `decostand' , NOTE: method = "pa"
otu_decostand <- decostand(randomized_otu, method = "pa")
# check total abundance in each sample
apply(otu_decostand, 1, sum)

## MBREJ515 MBREJ715 MBREJ915 MBREK515 MBREK715 MBREK915 MDPEJ515 MDPEJ715 MDPEJ915 MDPEK515 MDPEK715 MDPEK915 MINEJ515 MINEJ715 MINEJ915 MINEK515 MINEK715 MINEK915 MOTEJ515 MOTEJ715 MOTEJ915 MOTEK515 MOTEK715 MOTEK915 
##      906      574      434      268      256      358      493      447      476      276      284      381      840      632      586      452      383      511      505      343      444      274      238      381

# check for mismatches/missing species between community data and phylo tree
randomized_matches <- match.phylo.comm(randomized_tree, otu_decostand)
# the resulting object is a list with $phy and $comm elements.  replace our
# original data with the sorted/matched data
phy_randomized_rm2 <- randomized_matches$phy
comm_randomized_rm2 <- randomized_matches$comm

# Calculate the phylogenetic distances
phy_dist_randomized_rm2 <- cophenetic(phy_randomized_rm2)

  
## Calculate SES_MPD
###################################### INDEPENDENT SWAP ############################################
# calculate standardized effect size mean pairwise distance (ses.mpd)
unweighted_sesMPD_indepswap_randomized <- ses.mpd(comm_randomized_rm2, phy_dist_randomized_rm2, null.model = "independentswap", 
                                     abundance.weighted = FALSE, runs = 999)

df <- unweighted_sesMPD_indepswap_randomized

df$names <- row.names(df)
df_2 <- make_metadata_norep(df) %>%
  mutate(fraction = fct_recode(fraction, Particle = "WholePart", Free = "WholeFree"))
  
  
# Is there a relationship between richness and Unweighted Mean Pairwise distance?
summary(lm(ntaxa ~ mpd.obs.z, data = df_2))

## 
## Call:
## lm(formula = ntaxa ~ mpd.obs.z, data = df_2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -214.61 -128.20  -13.56   56.63  465.39 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   448.65      35.29  12.713  1.3e-11 ***
## mpd.obs.z      75.21      95.93   0.784    0.441    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 172.8 on 22 degrees of freedom
## Multiple R-squared:  0.02718,    Adjusted R-squared:  -0.01704 
## F-statistic: 0.6146 on 1 and 22 DF,  p-value: 0.4414

summary(lm(ntaxa ~ mpd.obs.z, data = filter(df_2, fraction == "Particle")))

## 
## Call:
## lm(formula = ntaxa ~ mpd.obs.z, data = filter(df_2, fraction == 
##     "Particle"))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -208.89 -109.98  -51.43   39.10  346.59 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   557.41      50.97  10.935 6.96e-07 ***
## mpd.obs.z     -18.70     153.76  -0.122    0.906    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 175.3 on 10 degrees of freedom
## Multiple R-squared:  0.001477,   Adjusted R-squared:  -0.09838 
## F-statistic: 0.01479 on 1 and 10 DF,  p-value: 0.9056

summary(lm(ntaxa ~ mpd.obs.z, data = filter(df_2, fraction == "Free")))

## 
## Call:
## lm(formula = ntaxa ~ mpd.obs.z, data = filter(df_2, fraction == 
##     "Free"))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -108.71  -51.56  -17.28   57.16  160.49 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   343.12      25.00  13.727 8.17e-08 ***
## mpd.obs.z      68.11      62.34   1.093      0.3    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 85.34 on 10 degrees of freedom
## Multiple R-squared:  0.1066, Adjusted R-squared:  0.01729 
## F-statistic: 1.194 on 1 and 10 DF,  p-value: 0.3002

randomized_rich_plot <- ggplot(df_2, aes(y = ntaxa, x = mpd.obs.z, fill = fraction)) +
  geom_vline(xintercept = 0, linetype = "dashed", size = 1.5) +
  geom_point(size = 3.5, aes(shape = season)) + ylab("Randomized Richness") +
  xlab("Unweighted Phylogenetic Diversity") +
  scale_fill_manual(values = fraction_colors) +
  scale_shape_manual(values = season_shapes) + 
  scale_x_continuous(limits = c(-1,1)) + 
  theme(legend.title = element_blank(), legend.position = "none")


plot_grid(randomized_rich_plot, season_legend,
          rel_heights = c(1, 0.05), nrow = 2, ncol = 1)

Figure S10

PCA with Euclidian Distances

# Make the Supplemental Figure
par(oma(0.1, 0.1, 0.1, 0.1))

## Error in oma(0.1, 0.1, 0.1, 0.1): could not find function "oma"

biplot(pca_environ,
     xlab = paste("PC1", paste(round(summary(pca_environ)$cont$importance[2,1]*100, digits = 2), "%", sep = ""), sep = ": "),
     ylab = paste("PC2", paste(round(summary(pca_environ)$cont$importance[2,2]*100, digits = 2), "%", sep = ""), sep = ": "))

par(mar = c(5,5,2,5))
plot(summary(pca_environ)$cont$importance[2,]*100, 
     xlab = "PCA Axis", 
     ylab = "Variation Explained Per Axis",
     ylim = c(0, 105),
     col = "cornflowerblue",
     cex =2,
     pch = 16)
par(new = T)
plot(summary(pca_environ)$cont$importance[3,]*100,
       cex = 2,
       pch = 17,
       ylim = c(0, 105),
       col = "firebrick3",
       axes=F, 
       xlab=NA, 
       ylab=NA)
axis(side = 4)
mtext(side = 4, line = 3, "Total Accumulated Variation")
legend("right",
       legend=c("Per Axis", "Total"),
       pch=c(16, 17), col=c("cornflowerblue", "firebrick3"))

The first two axes of the PCA explain 69.271% of the variation in the environmental data!

Linear Models with Environmental Variables

#####  SINGLE REGRESSION with PC1
## Free living
comm_lm_PC1_free <- summary(lm(frac_bacprod ~ PC1, 
           data = filter(metadata_pca, fraction == "Free")))
## Particle-associated
comm_lm_PC1_part <- summary(lm(frac_bacprod ~ PC1, 
           data = filter(metadata_pca, fraction == "Particle")))
## PER CAPITA: Free living
percap_lm_PC1_free <- summary(lm(log10(fracprod_per_cell_noinf) ~ PC1, 
           data = filter(metadata_pca, fraction == "Free")))
## PER CAPITA: Particle-associated
percap_lm_PC1_part <- summary(lm(log10(fracprod_per_cell_noinf) ~ PC1, 
           data = filter(metadata_pca, fraction == "Particle")))

# Put the results all together 
pca_lm_row1 <- c("Free", "Per-Liter","PC1",
          round(comm_lm_PC1_free$adj.r.squared, digits = 3),
          round(comm_lm_PC1_free$coefficients[2,4], digits = 3))
pca_lm_row2 <- c("Particle", "Per-Liter","PC1",
          round(comm_lm_PC1_part$adj.r.squared, digits = 3),
          round(comm_lm_PC1_part$coefficients[2,4], digits = 3))
pca_lm_row3 <- c("Free", "Per-Capita","PC1",
          round(percap_lm_PC1_free$adj.r.squared, digits = 3),
          round(percap_lm_PC1_free$coefficients[2,4], digits = 3))
pca_lm_row4 <- c("Particle", "Per-Capita", "PC1",
          round(percap_lm_PC1_part$adj.r.squared, digits = 3),
          round(percap_lm_PC1_part$coefficients[2,4], digits = 3))

pca_lm_results_df <- 
  rbind(pca_lm_row1, pca_lm_row2, pca_lm_row3, pca_lm_row4)

colnames(pca_lm_results_df) <- c("fraction", "Predictor","Prod_measure", "Adjusted_R2", "p-value")

row.names(pca_lm_results_df) = NULL

#####  SINGLE REGRESSION with PC2
## Free living
comm_lm_PC2_free <- summary(lm(frac_bacprod ~ PC2, 
           data = filter(metadata_pca, fraction == "Free")))
## Particle-associated
comm_lm_PC2_part <- summary(lm(frac_bacprod ~ PC2, 
           data = filter(metadata_pca, fraction == "Particle")))
## PER CAPITA: Free living
percap_lm_PC2_free <- summary(lm(log10(fracprod_per_cell_noinf) ~ PC2, 
           data = filter(metadata_pca, fraction == "Free")))
## PER CAPITA: Particle-associated
percap_lm_PC2_part <- summary(lm(log10(fracprod_per_cell_noinf) ~ PC2, 
           data = filter(metadata_pca, fraction == "Particle")))

# Put the results all together 
pca_lm2_row1 <- c("Free", "Per-Liter","PC2",
          round(comm_lm_PC2_free$adj.r.squared, digits = 3),
          round(comm_lm_PC2_free$coefficients[2,4], digits = 3))
pca_lm2_row2 <- c("Particle", "Per-Liter","PC2",
          round(comm_lm_PC2_part$adj.r.squared, digits = 3),
          round(comm_lm_PC2_part$coefficients[2,4], digits = 3))
pca_lm2_row3 <- c("Free", "Per-Capita","PC2",
          round(percap_lm_PC2_free$adj.r.squared, digits = 3),
          round(percap_lm_PC2_free$coefficients[2,4], digits = 3))
pca_lm2_row4 <- c("Particle", "Per-Capita", "PC2",
          round(percap_lm_PC2_part$adj.r.squared, digits = 3),
          round(percap_lm_PC2_part$coefficients[2,4], digits = 3))

pca_lm2_results_df <- 
  rbind(pca_lm2_row1, pca_lm2_row2, pca_lm2_row3, pca_lm2_row4)

colnames(pca_lm2_results_df) <- c("fraction", "Predictor","Prod_measure", "Adjusted_R2", "p-value")

row.names(pca_lm2_results_df) = NULL

#####  SINGLE REGRESSION with PC3
## Free living
comm_lm_PC3_free <- summary(lm(frac_bacprod ~ PC3, 
           data = filter(metadata_pca, fraction == "Free")))
## Particle-associated
comm_lm_PC3_part <- summary(lm(frac_bacprod ~ PC3, 
           data = filter(metadata_pca, fraction == "Particle")))
## PER CAPITA: Free living
percap_lm_PC3_free <- summary(lm(log10(fracprod_per_cell_noinf) ~ PC3, 
           data = filter(metadata_pca, fraction == "Free")))
## PER CAPITA: Particle-associated
percap_lm_PC3_part <- summary(lm(log10(fracprod_per_cell_noinf) ~ PC3, 
           data = filter(metadata_pca, fraction == "Particle")))

# Put the results all together 
pca_lm_row1_PC3 <- c("Free", "Per-Liter","PC3",
          round(comm_lm_PC3_free$adj.r.squared, digits = 3),
          round(comm_lm_PC3_free$coefficients[2,4], digits = 3))
pca_lm_row2_PC3 <- c("Particle", "Per-Liter","PC3",
          round(comm_lm_PC3_part$adj.r.squared, digits = 3),
          round(comm_lm_PC3_part$coefficients[2,4], digits = 3))
pca_lm_row3_PC3 <- c("Free", "Per-Capita","PC3",
          round(percap_lm_PC3_free$adj.r.squared, digits = 3),
          round(percap_lm_PC3_free$coefficients[2,4], digits = 3))
pca_lm_row4_PC3 <- c("Particle", "Per-Capita", "PC3",
          round(percap_lm_PC3_part$adj.r.squared, digits = 3),
          round(percap_lm_PC3_part$coefficients[2,4], digits = 3))

pca_lm_results_df_PC3 <- 
  rbind(pca_lm_row1_PC3, pca_lm_row2_PC3, pca_lm_row3_PC3, pca_lm_row4_PC3)

colnames(pca_lm_results_df_PC3) <- c("fraction", "Predictor","Prod_measure", "Adjusted_R2", "p-value")

row.names(pca_lm_results_df_PC3) = NULL

#####  SINGLE REGRESSION with PC4
## Free living
comm_lm_PC4_free <- summary(lm(frac_bacprod ~ PC4, 
           data = filter(metadata_pca, fraction == "Free")))
## Particle-associated
comm_lm_PC4_part <- summary(lm(frac_bacprod ~ PC4, 
           data = filter(metadata_pca, fraction == "Particle")))
## PER CAPITA: Free living
percap_lm_PC4_free <- summary(lm(log10(fracprod_per_cell_noinf) ~ PC4, 
           data = filter(metadata_pca, fraction == "Free")))
## PER CAPITA: Particle-associated
percap_lm_PC4_part <- summary(lm(log10(fracprod_per_cell_noinf) ~ PC4, 
           data = filter(metadata_pca, fraction == "Particle")))

# Put the results all together 
pca_lm_row1_PC4 <- c("Free", "Per-Liter","PC4",
          round(comm_lm_PC4_free$adj.r.squared, digits = 3),
          round(comm_lm_PC4_free$coefficients[2,4], digits = 3))
pca_lm_row2_PC4 <- c("Particle", "Per-Liter","PC4",
          round(comm_lm_PC4_part$adj.r.squared, digits = 3),
          round(comm_lm_PC4_part$coefficients[2,4], digits = 3))
pca_lm_row3_PC4 <- c("Free", "Per-Capita","PC4",
          round(percap_lm_PC4_free$adj.r.squared, digits = 3),
          round(percap_lm_PC4_free$coefficients[2,4], digits = 3))
pca_lm_row4_PC4 <- c("Particle", "Per-Capita", "PC4",
          round(percap_lm_PC4_part$adj.r.squared, digits = 3),
          round(percap_lm_PC4_part$coefficients[2,4], digits = 3))

pca_lm_results_df_PC4 <- 
  rbind(pca_lm_row1_PC4, pca_lm_row2_PC4, pca_lm_row3_PC4, pca_lm_row4_PC4)

colnames(pca_lm_results_df_PC4) <- c("fraction", "Predictor","Prod_measure", "Adjusted_R2", "p-value")

row.names(pca_lm_results_df_PC4) = NULL


#####  SINGLE REGRESSION with PC5
## Free living
comm_lm_PC5_free <- summary(lm(frac_bacprod ~ PC5, 
           data = filter(metadata_pca, fraction == "Free")))
## Particle-associated
comm_lm_PC5_part <- summary(lm(frac_bacprod ~ PC5, 
           data = filter(metadata_pca, fraction == "Particle")))
## PER CAPITA: Free living
percap_lm_PC5_free <- summary(lm(log10(fracprod_per_cell_noinf) ~ PC5, 
           data = filter(metadata_pca, fraction == "Free")))
## PER CAPITA: Particle-associated
percap_lm_PC5_part <- summary(lm(log10(fracprod_per_cell_noinf) ~ PC5, 
           data = filter(metadata_pca, fraction == "Particle")))

# Put the results all together 
pca_lm_row1_PC5 <- c("Free", "Per-Liter","PC5",
          round(comm_lm_PC5_free$adj.r.squared, digits = 3),
          round(comm_lm_PC5_free$coefficients[2,4], digits = 3))
pca_lm_row2_PC5 <- c("Particle", "Per-Liter","PC5",
          round(comm_lm_PC5_part$adj.r.squared, digits = 3),
          round(comm_lm_PC5_part$coefficients[2,4], digits = 3))
pca_lm_row3_PC5 <- c("Free", "Per-Capita","PC5",
          round(percap_lm_PC5_free$adj.r.squared, digits = 3),
          round(percap_lm_PC5_free$coefficients[2,4], digits = 3))
pca_lm_row4_PC5 <- c("Particle", "Per-Capita", "PC5",
          round(percap_lm_PC5_part$adj.r.squared, digits = 3),
          round(percap_lm_PC5_part$coefficients[2,4], digits = 3))

pca_lm_results_df_PC5 <- 
  rbind(pca_lm_row1_PC5, pca_lm_row2_PC5, pca_lm_row3_PC5, pca_lm_row4_PC5)

colnames(pca_lm_results_df_PC5) <- c("fraction", "Predictor","Prod_measure", "Adjusted_R2", "p-value")

row.names(pca_lm_results_df_PC5) = NULL


#####  SINGLE REGRESSION with PC6
## Free living
comm_lm_PC6_free <- summary(lm(frac_bacprod ~ PC6, 
           data = filter(metadata_pca, fraction == "Free")))
## Particle-associated
comm_lm_PC6_part <- summary(lm(frac_bacprod ~ PC6, 
           data = filter(metadata_pca, fraction == "Particle")))
## PER CAPITA: Free living
percap_lm_PC6_free <- summary(lm(log10(fracprod_per_cell_noinf) ~ PC6, 
           data = filter(metadata_pca, fraction == "Free")))
## PER CAPITA: Particle-associated
percap_lm_PC6_part <- summary(lm(log10(fracprod_per_cell_noinf) ~ PC6, 
           data = filter(metadata_pca, fraction == "Particle")))

# Put the results all together 
pca_lm_row1_PC6 <- c("Free", "Per-Liter","PC6",
          round(comm_lm_PC6_free$adj.r.squared, digits = 3),
          round(comm_lm_PC6_free$coefficients[2,4], digits = 3))
pca_lm_row2_PC6 <- c("Particle", "Per-Liter","PC6",
          round(comm_lm_PC6_part$adj.r.squared, digits = 3),
          round(comm_lm_PC6_part$coefficients[2,4], digits = 3))
pca_lm_row3_PC6 <- c("Free", "Per-Capita","PC6",
          round(percap_lm_PC6_free$adj.r.squared, digits = 3),
          round(percap_lm_PC6_free$coefficients[2,4], digits = 3))
pca_lm_row4_PC6 <- c("Particle", "Per-Capita", "PC6",
          round(percap_lm_PC6_part$adj.r.squared, digits = 3),
          round(percap_lm_PC6_part$coefficients[2,4], digits = 3))

pca_lm_results_df_PC6 <- 
  rbind(pca_lm_row1_PC6, pca_lm_row2_PC6, pca_lm_row3_PC6, pca_lm_row4_PC6)

colnames(pca_lm_results_df_PC6) <- c("fraction", "Predictor","Prod_measure", "Adjusted_R2", "p-value")

row.names(pca_lm_results_df_PC6) = NULL


################ Per-Capita production correlation with PC2
## 1. Run the linear regression 
lm_PC2_percap_part <- lm(log10(fracprod_per_cell_noinf) ~ PC2, data = dplyr::filter(metadata_pca, fraction == "Particle"))

# Linear models with PC1 and free-living PC2 are NOT SIGNIFICANT
summary(lm(log10(fracprod_per_cell_noinf) ~ PC2, data = dplyr::filter(metadata_pca, fraction == "Free")))

## 
## Call:
## lm(formula = log10(fracprod_per_cell_noinf) ~ PC2, data = dplyr::filter(metadata_pca, 
##     fraction == "Free"))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.6405 -0.1022 -0.0179  0.1564  0.5977 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -7.5750     0.1065 -71.098 7.39e-15 ***
## PC2          -0.1679     0.1048  -1.603     0.14    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3691 on 10 degrees of freedom
## Multiple R-squared:  0.2044, Adjusted R-squared:  0.1248 
## F-statistic: 2.569 on 1 and 10 DF,  p-value: 0.14

summary(lm(log10(fracprod_per_cell_noinf) ~ PC1, data = dplyr::filter(metadata_pca, fraction == "Free")))

## 
## Call:
## lm(formula = log10(fracprod_per_cell_noinf) ~ PC1, data = dplyr::filter(metadata_pca, 
##     fraction == "Free"))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.51287 -0.31915  0.00081  0.17747  0.65242 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -7.5750     0.1131 -67.006 1.34e-14 ***
## PC1           0.1199     0.1112   1.079    0.306    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3916 on 10 degrees of freedom
## Multiple R-squared:  0.1043, Adjusted R-squared:  0.0147 
## F-statistic: 1.164 on 1 and 10 DF,  p-value: 0.3059

summary(lm(log10(fracprod_per_cell_noinf) ~ PC1, data = dplyr::filter(metadata_pca, fraction == "Particle")))

## 
## Call:
## lm(formula = log10(fracprod_per_cell_noinf) ~ PC1, data = dplyr::filter(metadata_pca, 
##     fraction == "Particle"))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.64352 -0.35634 -0.01985  0.27133  1.09984 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -6.7108     0.1615 -41.563 1.35e-11 ***
## PC1           0.1568     0.1592   0.984    0.351    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5331 on 9 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.09719,    Adjusted R-squared:  -0.003121 
## F-statistic: 0.9689 on 1 and 9 DF,  p-value: 0.3507

# 2. Extract the R2 and p-value from the linear model: 
lm_lab_PC2_percap_part <- paste("atop(R^2 ==", round(summary(lm_PC2_percap_part)$adj.r.squared, digits = 3), ",",
            "p ==", round(unname(summary(lm_PC2_percap_part)$coefficients[,4][2]), digits = 3), ")")

## 3. Plot Per-Capita production correlation between particle and free
PC1_lm_plot <- 
  ggplot(metadata_pca, aes(x = PC1, y = log10(fracprod_per_cell_noinf), 
                         color = fraction,  fill = fraction)) +
  geom_point() + ylab("log10(Per-Capita Production)") + 
  scale_color_manual(values = fraction_colors) + 
  scale_fill_manual(values = fraction_colors) + 
  theme(legend.title = element_blank(), legend.position = "bottom")


PC2_lm_plot <- 
  ggplot(metadata_pca, aes(x = PC2, y = log10(fracprod_per_cell_noinf))) +
  geom_point(aes(color = fraction,  fill = fraction)) + ylab("log10(Per-Capita Production)") + 
  scale_color_manual(values = fraction_colors) + 
  scale_fill_manual(values = fraction_colors) + 
  geom_smooth(method = "lm", color = "#FF6600", fill = "#FF6600", data = dplyr::filter(metadata_pca, fraction == "Particle")) +
  # Add the R2 and p-value to the plot 
  annotate("text", y = -5.8, x=0.5, label=lm_lab_PC2_percap_part, parse = TRUE, color = "#FF6600", size = 4) +
  theme(legend.title = element_blank(), legend.position = "bottom")  

plot_grid(PC1_lm_plot, PC2_lm_plot, 
          labels = c("A", "B"),
          nrow = 1, ncol = 2)

Multiple Regression with PCA results

## Community-wide: Free living with PC1 + PC2 = Not significant
summary(lm(frac_bacprod ~ PC1 + PC2, data = filter(metadata_pca, fraction == "Free")))

## 
## Call:
## lm(formula = frac_bacprod ~ PC1 + PC2, data = filter(metadata_pca, 
##     fraction == "Free"))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -17.054  -9.327  -7.554  10.981  29.901 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   24.058      4.954   4.857   0.0009 ***
## PC1            6.295      4.871   1.292   0.2284    
## PC2           -4.385      4.871  -0.900   0.3914    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 17.16 on 9 degrees of freedom
## Multiple R-squared:  0.2161, Adjusted R-squared:  0.04186 
## F-statistic:  1.24 on 2 and 9 DF,  p-value: 0.3344

## Community-wide:Particle-associated with PC1 + PC2 = Not significant
summary(lm(frac_bacprod ~ PC1 + PC2, data = filter(metadata_pca, fraction == "Particle")))

## 
## Call:
## lm(formula = frac_bacprod ~ PC1 + PC2, data = filter(metadata_pca, 
##     fraction == "Particle"))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -12.481  -3.757   0.410   3.213  15.482 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)    9.958      2.129   4.678  0.00115 **
## PC1            1.107      2.093   0.529  0.60954   
## PC2           -4.424      2.093  -2.114  0.06368 . 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.373 on 9 degrees of freedom
## Multiple R-squared:  0.3454, Adjusted R-squared:  0.1999 
## F-statistic: 2.374 on 2 and 9 DF,  p-value: 0.1486

## PER CAPITA: Free living with PC1 + PC2 = Not significant
summary(lm(log10(fracprod_per_cell_noinf) ~ PC1 + PC2, data = filter(metadata_pca, fraction == "Free")))

## 
## Call:
## lm(formula = log10(fracprod_per_cell_noinf) ~ PC1 + PC2, data = filter(metadata_pca, 
##     fraction == "Free"))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.5064 -0.2177 -0.1287  0.2483  0.4932 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -7.5750     0.1047 -72.357  9.3e-14 ***
## PC1           0.1199     0.1029   1.165    0.274    
## PC2          -0.1679     0.1029  -1.631    0.137    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3627 on 9 degrees of freedom
## Multiple R-squared:  0.3087, Adjusted R-squared:  0.1551 
## F-statistic: 2.009 on 2 and 9 DF,  p-value: 0.1899

## PER CAPITA: Particle-associatedwith PC1 + PC2  = SIGNIFICANT (p = 0.015)
lm_percap_PC12_part <- lm(log10(fracprod_per_cell_noinf) ~ PC1 + PC2, data = filter(metadata_pca, fraction == "Particle"))
summary(lm_percap_PC12_part)

## 
## Call:
## lm(formula = log10(fracprod_per_cell_noinf) ~ PC1 + PC2, data = filter(metadata_pca, 
##     fraction == "Particle"))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.53476 -0.15558 -0.00926  0.08264  0.65269 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -6.6657     0.1074 -62.056 5.05e-12 ***
## PC1           0.2026     0.1060   1.912  0.09229 .  
## PC2          -0.3818     0.1075  -3.553  0.00747 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3522 on 8 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.6498, Adjusted R-squared:  0.5623 
## F-statistic: 7.423 on 2 and 8 DF,  p-value: 0.01504

Important! In a multiple regression, PC1 and PC2 correlate with log₁₀(per-capita heterotrophic production)! With an adjusted R² of 0.562 and a p-value of 0.015

Filtered PCA Linear Model Results

Below represents the significant LM results from the individual PCA axes:

# Put all of the PCA and environmental  linear model results together into one dataframe 
all_pca_lm_results <- dplyr::bind_rows(
  data.frame(pca_lm_results_df), 
  data.frame(pca_lm2_results_df), 
  data.frame(pca_lm_results_df_PC3),
  data.frame(pca_lm_results_df_PC4), 
  data.frame(pca_lm_results_df_PC5), 
  data.frame(pca_lm_results_df_PC6))

# Output the table in fancy form 
datatable(all_pca_lm_results, options = list(pageLength = 10)) 

# Output the results for models that have a pvalue of less than 0.05
pander(dplyr::filter(all_pca_lm_results, p.value < 0.05),
       caption = "Individual Linear Model results from the individual PCA axes with a p-value less than 0.06.")

Individual Linear Model results from the individual PCA axes with a p-value less than 0.06.

fraction	Predictor	Prod_measure	Adjusted_R2	p.value
Particle	Per-Capita	PC2	0.433	0.016

# Multiple regression with community-wide production
summary(lm(frac_bacprod ~ PC1 + PC2, data = filter(metadata_pca, fraction == "Particle")))

## 
## Call:
## lm(formula = frac_bacprod ~ PC1 + PC2, data = filter(metadata_pca, 
##     fraction == "Particle"))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -12.481  -3.757   0.410   3.213  15.482 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)    9.958      2.129   4.678  0.00115 **
## PC1            1.107      2.093   0.529  0.60954   
## PC2           -4.424      2.093  -2.114  0.06368 . 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.373 on 9 degrees of freedom
## Multiple R-squared:  0.3454, Adjusted R-squared:  0.1999 
## F-statistic: 2.374 on 2 and 9 DF,  p-value: 0.1486

summary(lm(Temp_C ~ PC2, data = filter(metadata_pca, fraction == "Particle")))

## 
## Call:
## lm(formula = Temp_C ~ PC2, data = filter(metadata_pca, fraction == 
##     "Particle"))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.0479 -1.2241 -0.8952  0.5037  4.1894 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  19.0575     0.6169  30.892 2.96e-11 ***
## PC2           3.0793     0.6066   5.076 0.000481 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.137 on 10 degrees of freedom
## Multiple R-squared:  0.7204, Adjusted R-squared:  0.6924 
## F-statistic: 25.77 on 1 and 10 DF,  p-value: 0.0004806

summary(lm(pH ~ PC2, data = filter(metadata_pca, fraction == "Particle")))

## 
## Call:
## lm(formula = pH ~ PC2, data = filter(metadata_pca, fraction == 
##     "Particle"))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.15389 -0.04228 -0.02834  0.03956  0.18650 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  8.31000    0.03184 261.003   <2e-16 ***
## PC2          0.09082    0.03131   2.901   0.0158 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1103 on 10 degrees of freedom
## Multiple R-squared:  0.4569, Adjusted R-squared:  0.4026 
## F-statistic: 8.414 on 1 and 10 DF,  p-value: 0.01581

ggplot(metadata_pca, aes(PC2, Temp_C)) + 
  geom_point() +
  geom_smooth(method = "lm")

ggplot(metadata_pca, aes(PC2, pH)) + 
  geom_point() +
  geom_smooth(method = "lm")

# Multiple regression with per-capita production
summary(lm(log10(fracprod_per_cell_noinf) ~ PC1 + PC2, data = filter(metadata_pca, fraction == "Free")))

## 
## Call:
## lm(formula = log10(fracprod_per_cell_noinf) ~ PC1 + PC2, data = filter(metadata_pca, 
##     fraction == "Free"))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.5064 -0.2177 -0.1287  0.2483  0.4932 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -7.5750     0.1047 -72.357  9.3e-14 ***
## PC1           0.1199     0.1029   1.165    0.274    
## PC2          -0.1679     0.1029  -1.631    0.137    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3627 on 9 degrees of freedom
## Multiple R-squared:  0.3087, Adjusted R-squared:  0.1551 
## F-statistic: 2.009 on 2 and 9 DF,  p-value: 0.1899

## Community Bulk Production: No correlation with PC2
summary(lm(frac_bacprod ~ Inverse_Simpson, data = lasso_data_df_particle)) 

## 
## Call:
## lm(formula = frac_bacprod ~ Inverse_Simpson, data = lasso_data_df_particle)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.5589 -2.1093 -0.1969  0.8752  7.8282 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     -0.45199    2.43866  -0.185  0.85666    
## Inverse_Simpson  0.29344    0.05781   5.076  0.00048 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.571 on 10 degrees of freedom
## Multiple R-squared:  0.7204, Adjusted R-squared:  0.6925 
## F-statistic: 25.77 on 1 and 10 DF,  p-value: 0.0004804

summary(lm(frac_bacprod ~ PC2, data = lasso_data_df_particle)) 

## 
## Call:
## lm(formula = frac_bacprod ~ PC2, data = lasso_data_df_particle)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -11.5929  -4.0916   0.8315   2.6669  16.4468 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    9.958      2.050   4.857 0.000665 ***
## PC2           -4.424      2.016  -2.194 0.052954 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.103 on 10 degrees of freedom
## Multiple R-squared:  0.325,  Adjusted R-squared:  0.2575 
## F-statistic: 4.815 on 1 and 10 DF,  p-value: 0.05295

## Per Capita Production --> Correlations for both Simpson and PC2
summary(lm(log10(fracprod_per_cell_noinf) ~ Inverse_Simpson, data = lasso_data_df_particle)) 

## 
## Call:
## lm(formula = log10(fracprod_per_cell_noinf) ~ Inverse_Simpson, 
##     data = lasso_data_df_particle)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.28651 -0.18384 -0.11125  0.07337  0.56444 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     -7.360961   0.159490 -46.153 5.27e-12 ***
## Inverse_Simpson  0.018087   0.003759   4.812 0.000958 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2969 on 9 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.7201, Adjusted R-squared:  0.689 
## F-statistic: 23.15 on 1 and 9 DF,  p-value: 0.0009581

summary(lm(log10(fracprod_per_cell_noinf) ~ PC2, data = lasso_data_df_particle)) 

## 
## Call:
## lm(formula = log10(fracprod_per_cell_noinf) ~ PC2, data = lasso_data_df_particle)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.48197 -0.30255  0.04728  0.12955  0.87488 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -6.6877     0.1215  -55.03 1.09e-12 ***
## PC2          -0.3568     0.1214   -2.94   0.0165 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4008 on 9 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.4898, Adjusted R-squared:  0.4332 
## F-statistic: 8.642 on 1 and 9 DF,  p-value: 0.01649

summary(lm(log10(fracprod_per_cell_noinf) ~ Inverse_Simpson + PC2, data = lasso_data_df_particle)) 

## 
## Call:
## lm(formula = log10(fracprod_per_cell_noinf) ~ Inverse_Simpson + 
##     PC2, data = lasso_data_df_particle)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.33410 -0.17773 -0.06805  0.05367  0.53447 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     -7.264550   0.234527 -30.975 1.28e-09 ***
## Inverse_Simpson  0.015588   0.005811   2.682   0.0278 *  
## PC2             -0.080742   0.138992  -0.581   0.5773    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3084 on 8 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.7314, Adjusted R-squared:  0.6643 
## F-statistic: 10.89 on 2 and 8 DF,  p-value: 0.005205

summary(lm(log10(fracprod_per_cell_noinf) ~ unweighted_div, data = lasso_data_df_particle)) 

## Error in eval(predvars, data, env): object 'unweighted_div' not found

Environmental Params: as they are

metadata_pca_particle <- metadata_pca %>%
  filter(fraction == "Particle")

# Temperature
summary(lm(frac_bacprod ~ Temp_C, data = metadata_pca_particle))

## 
## Call:
## lm(formula = frac_bacprod ~ Temp_C, data = metadata_pca_particle)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -11.775  -4.172  -1.356   3.410  15.237 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  30.6913    11.3062   2.715   0.0218 *
## Temp_C       -1.0879     0.5825  -1.868   0.0913 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.444 on 10 degrees of freedom
## Multiple R-squared:  0.2586, Adjusted R-squared:  0.1845 
## F-statistic: 3.489 on 1 and 10 DF,  p-value: 0.09133

summary(lm(log10(fracprod_per_cell_noinf) ~ Temp_C, data = metadata_pca_particle))

## 
## Call:
## lm(formula = log10(fracprod_per_cell_noinf) ~ Temp_C, data = metadata_pca_particle)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.46741 -0.21600  0.05621  0.11958  0.67225 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -4.63630    0.62565  -7.410 4.06e-05 ***
## Temp_C      -0.10720    0.03158  -3.395  0.00794 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3716 on 9 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.5615, Adjusted R-squared:  0.5128 
## F-statistic: 11.52 on 1 and 9 DF,  p-value: 0.007941

summary(lm(log10(fracprod_per_cell_noinf) ~ Temp_C, data = metadata_pca_particle))

## 
## Call:
## lm(formula = log10(fracprod_per_cell_noinf) ~ Temp_C, data = metadata_pca_particle)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.46741 -0.21600  0.05621  0.11958  0.67225 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -4.63630    0.62565  -7.410 4.06e-05 ***
## Temp_C      -0.10720    0.03158  -3.395  0.00794 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3716 on 9 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.5615, Adjusted R-squared:  0.5128 
## F-statistic: 11.52 on 1 and 9 DF,  p-value: 0.007941

# pH
summary(lm(frac_bacprod ~ pH, data = metadata_pca_particle))

## 
## Call:
## lm(formula = frac_bacprod ~ pH, data = metadata_pca_particle)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -10.0815  -3.2449  -0.2237   2.3649   9.9338 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)   380.72      96.43   3.948  0.00274 **
## pH            -44.62      11.60  -3.845  0.00324 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.491 on 10 degrees of freedom
## Multiple R-squared:  0.5966, Adjusted R-squared:  0.5562 
## F-statistic: 14.79 on 1 and 10 DF,  p-value: 0.003237

summary(lm(log10(fracprod_per_cell_noinf) ~ pH, data = metadata_pca_particle))

## 
## Call:
## lm(formula = log10(fracprod_per_cell_noinf) ~ pH, data = metadata_pca_particle)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.3969 -0.2459 -0.1015  0.2103  0.5372 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)   17.126      6.249   2.741  0.02283 * 
## pH            -2.867      0.751  -3.817  0.00411 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3467 on 9 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.6182, Adjusted R-squared:  0.5758 
## F-statistic: 14.57 on 1 and 9 DF,  p-value: 0.004106

# SO4_mgperL
summary(lm(frac_bacprod ~ SO4_mgperL, data = metadata_pca_particle))

## 
## Call:
## lm(formula = frac_bacprod ~ SO4_mgperL, data = metadata_pca_particle)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -7.477 -5.981 -1.670  2.509 20.736 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  8.55695   13.79530   0.620    0.549
## SO4_mgperL   0.08713    0.84360   0.103    0.920
## 
## Residual standard error: 8.641 on 10 degrees of freedom
## Multiple R-squared:  0.001066,   Adjusted R-squared:  -0.09883 
## F-statistic: 0.01067 on 1 and 10 DF,  p-value: 0.9198

summary(lm(log10(fracprod_per_cell_noinf) ~ SO4_mgperL, data = metadata_pca_particle))

## 
## Call:
## lm(formula = log10(fracprod_per_cell_noinf) ~ SO4_mgperL, data = metadata_pca_particle)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.56552 -0.38872 -0.06665  0.13798  1.22608 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -6.38030    0.90018  -7.088 5.74e-05 ***
## SO4_mgperL  -0.02135    0.05465  -0.391    0.705    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5564 on 9 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.01667,    Adjusted R-squared:  -0.09258 
## F-statistic: 0.1526 on 1 and 9 DF,  p-value: 0.7051

Bonus Visualizations/Analysis

Community-wide and Per-Capita Heterotrophic Production over Station and Season

# Plot the community wide heterotrophic productivity over station and season
prod1 <- ggplot(dplyr::filter(metadata, year == "2015" & fraction %in% c("Free", "Particle")), 
       aes(x = lakesite, y = frac_bacprod, fill = fraction)) + 
  geom_bar(stat = "identity", color = "black",  position=position_dodge()) + 
  geom_errorbar(aes(ymin = frac_bacprod - SD_frac_bacprod, ymax = frac_bacprod + SD_frac_bacprod), 
                position = position_dodge(width = 0.75), width = 0.25) +
  facet_grid(. ~ season) + 
  scale_y_continuous(expand = c(0,0), limits = c(0, 80), 
                     breaks = seq(from = 0, to = 80, by = 15)) +
  scale_fill_manual(values = fraction_colors) + 
  guides(fill=guide_legend(nrow=1,byrow=TRUE)) +
  xlab("Lake Station") + ylab("Community Production\n(ug C/L/day)") + 
  theme(axis.text.x = element_blank(), #element_text(angle = 30, vjust = 1, hjust = 1),
        axis.title.x = element_blank(),
        legend.position = c(0.65, 0.9), legend.title = element_blank(),
        plot.margin = unit(c(0.25, 0.25, 0.25, 0.25), "cm")) #top, right, bottom, and left margins)

# Plot the per-capita heterotrophic productivity over station and season
prod2 <- ggplot(dplyr::filter(metadata, year == "2015" & fraction %in% c("Free", "Particle")), 
       aes(x = lakesite, y = fracprod_per_cell_noinf, fill = fraction)) + 
  geom_bar(stat = "identity", color = "black",  position=position_dodge()) + 
  facet_grid(. ~ season) +
  scale_y_continuous(expand = c(0,0)) +
  scale_fill_manual(values = fraction_colors) + 
  guides(fill=guide_legend(nrow=1,byrow=TRUE)) +
  xlab("Lake Station") + ylab("Per-Capita Production \n(ug C/cell/day)") + 
  theme(axis.text.x = element_text(angle = 30, vjust = 1, hjust = 1),
        legend.position = c(0.65, 0.9), legend.title = element_blank(),
        plot.margin = unit(c(0.25, 0.25, 0.25, 0.5), "cm")) #top, right, bottom, and left margins)

# Put the two plots together 
plot_grid(prod1, prod2, labels = c("A", "B"),
          rel_heights = c(1, 1.3),
          nrow = 2, ncol = 1,
          align = "v")

Total Production

### Total production by station 
prod_by_station <- ggplot(dplyr::filter(metadata, year == "2015" & fraction %in% c("Free")), 
       aes(x = lakesite, y = tot_bacprod)) + 
  geom_boxplot(fill = "grey") + geom_point(size = 3.5, position = position_jitter()) +
  scale_y_continuous(expand = c(0,0), limits = c(0,100)) +
  scale_fill_manual(values = fraction_colors) + 
  guides(fill=guide_legend(nrow=1,byrow=TRUE)) +
  xlab("Lake Station") + ylab("Total Production (ug C/cell/day)") + 
  theme(axis.text.x = element_text(angle = 30, vjust = 1, hjust = 1),
        legend.position = c(0.85, 0.9), legend.title = element_blank(),
        plot.margin = unit(c(0.25, 0.25, 0.25, 0.5), "cm")) #top, right, bottom, and left margins)

### Total production by Season 
prod_by_season <- ggplot(dplyr::filter(metadata, year == "2015" & fraction %in% c("Free")), 
       aes(x = season, y = tot_bacprod)) + 
  geom_boxplot(fill = "grey") + geom_point(size = 3.5, position = position_jitter()) +
  scale_y_continuous(expand = c(0,0), limits = c(0,100)) +
  scale_fill_manual(values = fraction_colors) + 
  guides(fill=guide_legend(nrow=1,byrow=TRUE)) +
  xlab("Season") + 
  theme(axis.text.x = element_text(angle = 30, vjust = 1, hjust = 1),
        axis.title.y = element_blank(),
        legend.position = c(0.85, 0.9), legend.title = element_blank(),
        plot.margin = unit(c(0.25, 0.25, 0.25, 0.5), "cm")) #top, right, bottom, and left margins)

plot_grid(prod_by_station, prod_by_season, labels = c("A", "B"),
          rel_widths = c(1, 0.75),
          nrow = 1, ncol = 2,
          align = "h")

The highest total heterotrophic production rates occur in the river station and in the spring and summer.

Environmental Variables Pairs

pairs(environmental_data)

Session Information

devtools::session_info() # This will include session info with all R package version information

##  setting  value                       
##  version  R version 3.4.2 (2017-09-28)
##  system   x86_64, darwin15.6.0        
##  ui       X11                         
##  language (EN)                        
##  collate  en_US.UTF-8                 
##  tz       America/Detroit             
##  date     2017-11-22                  
## 
##  package      * version  date       source        
##  ade4           1.7-8    2017-08-09 CRAN (R 3.4.1)
##  ape          * 4.1      2017-02-14 CRAN (R 3.4.0)
##  assertthat     0.2.0    2017-04-11 CRAN (R 3.4.0)
##  backports      1.1.1    2017-09-25 CRAN (R 3.4.2)
##  base         * 3.4.2    2017-10-04 local         
##  bindr          0.1      2016-11-13 CRAN (R 3.4.0)
##  bindrcpp     * 0.2      2017-06-17 CRAN (R 3.4.0)
##  Biobase        2.30.0   2017-10-05 Bioconductor  
##  BiocGenerics   0.22.1   2017-10-07 Bioconductor  
##  biomformat     1.4.0    2017-04-25 Bioconductor  
##  Biostrings     2.44.2   2017-07-21 Bioconductor  
##  broom        * 0.4.2    2017-02-13 CRAN (R 3.4.0)
##  car          * 2.1-5    2017-07-04 CRAN (R 3.4.1)
##  caret        * 6.0-77   2017-09-07 CRAN (R 3.4.1)
##  class          7.3-14   2015-08-30 CRAN (R 3.4.2)
##  cluster        2.0.6    2017-03-10 CRAN (R 3.4.2)
##  codetools      0.2-15   2016-10-05 CRAN (R 3.4.2)
##  colorspace     1.3-2    2016-12-14 CRAN (R 3.4.0)
##  compiler       3.4.2    2017-10-04 local         
##  cowplot      * 0.8.0    2017-07-30 CRAN (R 3.4.1)
##  CVST           0.2-1    2013-12-10 CRAN (R 3.4.0)
##  data.table     1.10.4   2017-02-01 CRAN (R 3.4.0)
##  datasets     * 3.4.2    2017-10-04 local         
##  ddalpha        1.3.1    2017-09-27 CRAN (R 3.4.2)
##  DEoptimR       1.0-8    2016-11-19 CRAN (R 3.4.0)
##  devtools     * 1.13.3   2017-08-02 CRAN (R 3.4.1)
##  digest         0.6.12   2017-01-27 CRAN (R 3.4.0)
##  dimRed         0.1.0    2017-05-04 CRAN (R 3.4.0)
##  dplyr        * 0.7.4    2017-09-28 CRAN (R 3.4.2)
##  DRR            0.0.2    2016-09-15 CRAN (R 3.4.0)
##  DT           * 0.2      2016-08-09 CRAN (R 3.4.0)
##  evaluate       0.10.1   2017-06-24 CRAN (R 3.4.1)
##  forcats      * 0.2.0    2017-01-23 CRAN (R 3.4.0)
##  foreach      * 1.4.3    2015-10-13 CRAN (R 3.4.0)
##  foreign        0.8-69   2017-06-22 CRAN (R 3.4.2)
##  geosphere      1.5-5    2016-06-15 CRAN (R 3.4.0)
##  ggmap        * 2.6.1    2016-01-23 CRAN (R 3.4.0)
##  ggplot2      * 2.2.1    2016-12-30 CRAN (R 3.4.0)
##  glmnet       * 2.0-13   2017-09-22 CRAN (R 3.4.2)
##  glue           1.1.1    2017-06-21 cran (@1.1.1) 
##  gower          0.1.2    2017-02-23 CRAN (R 3.4.0)
##  graphics     * 3.4.2    2017-10-04 local         
##  grDevices    * 3.4.2    2017-10-04 local         
##  grid           3.4.2    2017-10-04 local         
##  gtable         0.2.0    2016-02-26 CRAN (R 3.4.0)
##  htmltools      0.3.6    2017-04-28 CRAN (R 3.4.0)
##  htmlwidgets    0.9      2017-07-10 CRAN (R 3.4.1)
##  igraph         1.1.2    2017-07-21 CRAN (R 3.4.1)
##  ipred          0.9-6    2017-03-01 CRAN (R 3.4.0)
##  IRanges        2.10.5   2017-10-08 Bioconductor  
##  iterators      1.0.8    2015-10-13 CRAN (R 3.4.0)
##  jpeg           0.1-8    2014-01-23 CRAN (R 3.4.0)
##  jsonlite       1.5      2017-06-01 CRAN (R 3.4.0)
##  kernlab        0.9-25   2016-10-03 CRAN (R 3.4.0)
##  knitr          1.17     2017-08-10 CRAN (R 3.4.1)
##  labeling       0.3      2014-08-23 CRAN (R 3.4.0)
##  lattice      * 0.20-35  2017-03-25 CRAN (R 3.4.2)
##  lava           1.5.1    2017-09-27 CRAN (R 3.4.2)
##  lazyeval       0.2.0    2016-06-12 CRAN (R 3.4.0)
##  lme4           1.1-14   2017-09-27 CRAN (R 3.4.2)
##  lubridate      1.6.0    2016-09-13 CRAN (R 3.4.0)
##  magrittr       1.5      2014-11-22 CRAN (R 3.4.0)
##  mapproj        1.2-5    2017-06-08 CRAN (R 3.4.0)
##  maps           3.2.0    2017-06-08 CRAN (R 3.4.0)
##  MASS         * 7.3-47   2017-02-26 CRAN (R 3.4.2)
##  Matrix       * 1.2-11   2017-08-21 CRAN (R 3.4.2)
##  MatrixModels   0.4-1    2015-08-22 CRAN (R 3.4.0)
##  memoise        1.1.0    2017-04-21 CRAN (R 3.4.0)
##  methods      * 3.4.2    2017-10-04 local         
##  mgcv           1.8-22   2017-09-19 CRAN (R 3.4.2)
##  minqa          1.2.4    2014-10-09 CRAN (R 3.4.0)
##  mnormt         1.5-5    2016-10-15 CRAN (R 3.4.0)
##  ModelMetrics   1.1.0    2016-08-26 CRAN (R 3.4.0)
##  multtest       2.32.0   2017-04-25 Bioconductor  
##  munsell        0.4.3    2016-02-13 CRAN (R 3.4.0)
##  nlme         * 3.1-131  2017-02-06 CRAN (R 3.4.2)
##  nloptr         1.0.4    2014-08-04 CRAN (R 3.4.0)
##  nnet           7.3-12   2016-02-02 CRAN (R 3.4.2)
##  pander       * 0.6.1    2017-08-06 CRAN (R 3.4.1)
##  parallel       3.4.2    2017-10-04 local         
##  pbkrtest       0.4-7    2017-03-15 CRAN (R 3.4.0)
##  permute      * 0.9-4    2016-09-09 CRAN (R 3.4.0)
##  phyloseq     * 1.20.0   2017-04-25 Bioconductor  
##  picante      * 1.6-2    2014-03-05 CRAN (R 3.4.0)
##  pkgconfig      2.0.1    2017-03-21 CRAN (R 3.4.0)
##  plyr           1.8.4    2016-06-08 CRAN (R 3.4.0)
##  png            0.1-7    2013-12-03 CRAN (R 3.4.0)
##  prodlim        1.6.1    2017-03-06 CRAN (R 3.4.0)
##  proto          1.0.0    2016-10-29 CRAN (R 3.4.0)
##  psych          1.7.8    2017-09-09 CRAN (R 3.4.2)
##  purrr          0.2.4    2017-10-18 CRAN (R 3.4.2)
##  quantreg       5.33     2017-04-18 CRAN (R 3.4.0)
##  R6             2.2.2    2017-06-17 CRAN (R 3.4.0)
##  Rcpp           0.12.13  2017-09-28 CRAN (R 3.4.2)
##  RcppRoll       0.2.2    2015-04-05 CRAN (R 3.4.0)
##  recipes        0.1.0    2017-07-27 CRAN (R 3.4.1)
##  reshape2       1.4.2    2016-10-22 CRAN (R 3.4.0)
##  RgoogleMaps    1.4.1    2016-09-18 CRAN (R 3.4.0)
##  rhdf5          2.20.0   2017-04-25 Bioconductor  
##  rjson          0.2.15   2014-11-03 CRAN (R 3.4.0)
##  rlang          0.1.2    2017-08-09 cran (@0.1.2) 
##  rmarkdown      1.6      2017-06-15 CRAN (R 3.4.0)
##  robustbase     0.92-7   2016-12-09 CRAN (R 3.4.0)
##  rpart          4.1-11   2017-03-13 CRAN (R 3.4.2)
##  rprojroot      1.2      2017-01-16 CRAN (R 3.4.0)
##  S4Vectors      0.14.7   2017-10-08 Bioconductor  
##  sandwich     * 2.4-0    2017-07-26 CRAN (R 3.4.1)
##  scales         0.5.0    2017-08-24 CRAN (R 3.4.1)
##  sfsmisc        1.1-1    2017-06-08 CRAN (R 3.4.0)
##  sp             1.2-5    2017-06-29 CRAN (R 3.4.1)
##  SparseM        1.77     2017-04-23 CRAN (R 3.4.0)
##  splines        3.4.2    2017-10-04 local         
##  stats        * 3.4.2    2017-10-04 local         
##  stats4         3.4.2    2017-10-04 local         
##  stringi        1.1.5    2017-04-07 CRAN (R 3.4.0)
##  stringr        1.2.0    2017-02-18 CRAN (R 3.4.0)
##  survival       2.41-3   2017-04-04 CRAN (R 3.4.2)
##  tibble       * 1.3.4    2017-08-22 CRAN (R 3.4.1)
##  tidyr        * 0.7.2    2017-10-16 CRAN (R 3.4.2)
##  tidyselect     0.2.2    2017-10-10 CRAN (R 3.4.2)
##  timeDate       3012.100 2015-01-23 CRAN (R 3.4.0)
##  tools          3.4.2    2017-10-04 local         
##  utils        * 3.4.2    2017-10-04 local         
##  vegan        * 2.4-4    2017-08-24 CRAN (R 3.4.1)
##  withr          2.0.0    2017-07-28 CRAN (R 3.4.1)
##  XVector        0.16.0   2017-04-25 Bioconductor  
##  yaml           2.1.14   2016-11-12 CRAN (R 3.4.0)
##  zlibbioc       1.16.0   2017-10-05 Bioconductor  
##  zoo            1.8-0    2017-04-12 CRAN (R 3.4.0)